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Summary  of  Work  Completed 

Research  performed  under  this  grant  contributed  to  papers  [1]— [8]  listed  below  in  the 

references.  The  papers  fall  into  four  areas. 

1.  [1]— [3]  concerned  arbitrarily  varying  channels.  Preliminary  versions  of  [1]  and  [3]  were 
included  in  last  year’s  report.  These  papers  have  since  appeared  in  the  literature  and 
are  therefore  not  included  in  this  report. 

2.  [4]  and  [8]  focused  on  the  design  of  distributed  estimation  systems  subject  to  communi¬ 
cation  and  computation  constraints.  The  theory  can  be  found  in  [4]  and  a  description 
and  listing  of  the  programs  used  to  implement  the  design  algorithm  can  be  found  in 
[8].  Both  of  these  papers  are  included  with  this  report. 

3.  [5]  and  [7]  concerned  the  computation  of  shot-noise  probability  distributions  and  image 
detection  based  on  shot-noise  observations.  These  papers  are  included  with  this  report. 

4.  [6]  is  a  tutorial  paper  on  wavelet  transforms  for  discrete-time  periodic  signals.  This 
paper  is  included  with  this  report. 
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Distributed  Estimation  and  Quantization 

John  A.  Gubner.  Member.  IEEE 


Abstract  —  We  develop  An  algorithm  for  the  design  of  an  n- 
sensor  distributed  estimation  system  subject  to  communication 
and  computation  constraints.  The  algorithm  uses  only  bivariate 
probability  distributions  and  yields  locally  optimal  estimators 
that  satisfy  the  required  system  constraints.  It  is  shown  that 
the  algorithm  is  a  generalization  of  the  classical  Lloyd-Max 
results. 

Index  Terms  —  Distributed  estimation,  quantization,  Lloyd- 
Max  algorithm. 

I.  Introduction 

Consider  the  distributed  estimation  system  shown  in 
Fig.  1.  The  system  consists  of  n  sensor  platforms  whose 
respective  measurements,  Vj,...,V'n,  are  related  to  some 
unobservable  quantity,  say  A".  Each  sensor  platform  pro¬ 
cesses  its  respective  measurement  and  transmits  the  result 
over  a  communication  channel  to  a  common  fusion  cen¬ 
ter.  The  sensors  do  not  communicate  with  each  other,  and 
there  is  no  feedback  from  the  fusion  center  to  the  sensor 
platforms.  The  task  of  the  fusion  center  is  to  estimate  the 
unobserved  quantity  X.  We  denote  this  estimate  by  ,Y. 
Clearly,  X  is  a  function  of  Yi , . . . ,  Y„ ,  and  we  can  write 

X  =  /(Yi . Y„)  for  some  function  /.  The  problem  then 

is  to  choose  the  function  /  so  that  X  is  close  to  X  in  some 
sense.  For  example,  it  is  well  known  that  in  the  appropri¬ 
ate  probabilistic  setting,  the  minimum-mean-square-error 
estimate  of  ,Y  given  Yi . Y„  is  the  conditional  expec¬ 

tation  of  ,Y  given  Yi .....  Y„,  denoted  B[X  |  Yj, . . Y„  ]. 
However,  there  are  many  situations  in  which  the  condi¬ 
tional  expectation  does  not  provide  a  satisfactory  solution 
to  the  problem  of  choosing  /: 

1 )  In  general,  the  functional  form  of  E[  X  |  Yi , . . . ,  Y„  ]  as 
a  function  of  Yi , . . . ,  Y„  is  difficult  determine,  and  it 
requires  knowledge  of  the  joint  probability  distribu¬ 
tion  of  ,Y,  Yi, . . . ,  Yn.  In  practice  this  complete  joint 
distribution  may  not  be  available. 

2)  To  compute  E[  X  (  Ylt . . . ,  Y„  ],  the  fusion  center  must 

in  general  have  access  to  all  of  the  sensor  measure¬ 
ments  Yi , . . . ,  Y„.  Hence,  even  if  the  sensor  platforms 
have  local  processing  capability,  it  is  of  little  use  in 
computing  E[,Y  |  Y\ . Y„],  If  the  number  of  sen¬ 

sor  platforms  is  very  large,  the  burden  of  computing 


This  research  wag  presented  in  part  at  the  1990  IEEE  International 
Symposium  on  Information  Theory,  San  Diego,  CA,  Jan.  14-19,  1990. 
This  work  wag  supported  in  part  by  the  Air  Force  Office  of  Scientific 
Research  under  Grant  AFOSR-90-0181. 

The  author  is  with  the  Department  of  Electrical  and  Computer 
Engineering,  University  of  Wisconsin,  Madison,  WI  53706. 


E[ X  (  Yi , ....  Y„  ]  at  the  fusion  center,  even  if  the  for¬ 
mula  is  relatively  simple,  may  be  prohibitive.  Such 
considerations  are  important  if  the  estimate  of  X  must 
be  computed  in  real  time.  By  using  a  suboptimal  es¬ 
timator  of  X  for  which  some  of  the  processing  can  be 
done  locally  at  the  sensor  platforms,  it  may  be  possi¬ 
ble  to  design  an  acceptable  estimator  that  can  operate 
in  real  time. 

3)  As  indicated  in  Fig.  1,  the  sensor  platforms  transmit 
their  data  to  the  fusion  center.  However,  using  any 
physical  communication  system,  it  is  not  possible  to 
transmit  real- valued  quantities  without  distortion.  In 
this  situation,  the  conditional  expectation,  or  even  the 
best  linear  estimate,  is  generally  a  physically  unreal¬ 
izable  solution 

In  this  paper  we  develop  an  algorithm  to  design  solutions 
to  the  distributed  estimation  problem  that  do  not  suffer 
from  these  difficulties. 


Vj  Y„ 


Estimate  of  X 


Fig.  1.  A  Distributed  Estimation  System. 


II.  Background  and  Notation 

Our  approach  is  to  consider  quantization  for  distributed 
estimation  systems.  The  goal  of  quantization  in  such  sys¬ 
tems  is  to  provide  a  good  estimate  of  the  unobservable, 
X,  rather  than  to  reconstruct  the  sensor  measurements 
Yi, . . . ,  Yn  as  in  [3].  Quantization  for  estimation  has  been 
studied  for  a  single  sensor  by  Ephraim  and  Gray  [2]  and 
by  Ayanoglu  (1],  The  multi-sensor  case  has  been  studied 
by  Lam  and  Reibman  [5],  and  we  discuss  their  work  in 
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more  detail  below.  The  paper  by  Zhang  and  Berger  [9] 
considers  an  asymptotic  estimation  problem  in  which  the 
observations  are  discrete  random  variables  taking  finitely 
many  values  and  the  unobservable  quantity  is  not  a  ran¬ 
dom  variable,  but  a  deterministic  and  unknown  parameter 
in  some  finite-dimensional  Euclidean  space. 

.4.  System  Model 

Let  ,Y.  Vi, . . . ,  Yn  be  real- valued  random  variables  on 
some  probability  space  (0,-F,  P).  Each  sensor  platform 
k  processes  its  measurement  Yt  to  obtain  an  output  2*. 
Each  Zt  is  then  transmitted  to  the  fusion  center.  We 
assume  that  the  communication  channel  connecting  the 
sensors  to  the  fusion  center  has  a  positive  capacity,  and 
that  the  use  of  error-correcting  codes  permits  us  to  view 
the  channel  as  noiseless.  We  suppose  that  the  channel 
ran  transmit  messages  of  log2  .V  bits  without  error,  where 
.V  >  2  is  an  integer.  For  each  k,  let  .4*1 . 4*.v  be  a  par¬ 

tition  of  IR.  We  require  that  the  sensor  platform  output 
Zk  be  given  by 


Clearly,  in  order  to  compute  (3),  we  need  to  know 
E[  X  |  Vi , . . . ,  Yn  ]•  If  the  entire  joint  distribution  Fxy,  y„ 
is  not  available,  computation  of  h  will  not  be  possible  in 
general.  Another  consideration  in  some  applications  is  the 
computation  of  (3)  in  real  time.  If  (3)  is  not  computable  in 
real  time,  all  the  different  possible  values  of  the  right-hand 
side  of  (3)  will  have  to  be  precomputed  and  stored.  For 
an  n-sensor  system  with  N-component  partitions,  there 
are  Nn  different  numbers  to  compute  and  store.  Finally, 
if  more  sensors  are  added  at  a  later  date,  there  will  be  no 
way  to  take  advantage  of  the  work  already  done  to  develop 
the  n-sensor  system;  all  of  the  numbers  given  by  (3)  will 
have  to  be  recomputed  for  an  even  larger  value  of  n. 

The  preceding  paragraph  assumed  that  the  partitions 
were  given.  If  h  is  arbitrary  and  given,  and  the  partitions 
are  given  for  k  ^  /,  then  the  remaining  partition 
should  satisfy  (in  order  to  minimize  (2)  [5]) 

y  £  Ai,  <=>  V,,(y)  <V,}{y)  for  all  j  =  1, - ,V 

(4) 

where, 


,v 

2*  =  nu.Ot).  (i) 

■  =  l 

where  lily)  denotes  the  indicator  function  of  the  set  .4  C 
IR;  i.e  .  l.\(y)  =  l  if  y  €  .4  and  l.t(y)  —  0  otherwise. 

Under  the  preceding  constraints,  the  function  /  dis¬ 
cussed  in  Section  I  must  be  of  the  form 

/(V'i . V„)  =  h(Z , . Zn), 

where  each  Z*  is  equal  to  the  function  of  Yt  determined 

by  ( 1 ). 

B.  Relation  to  [5] 

We  now  briefly  summarize  the  approach  in  (5].  If  the  sets 
(.4*, }  are  fixed,  one  wants  to  find  a  function  h(Z\ , . . . ,  Z„) 
that  minimizes  the  mean-square  error, 


U,.(y)  =  E[|E[X|V\,...,Y'n]-M‘-l)l2|Vj  =  y], 

and  /»/(»  -  1)  =  h(Zi . -  l,Z|+i,. .  .  ,Z„).  The 

approach  in  [5]  was  to  use  (3)  and  (4)  as  the  basis  of  an 
algorithm  for  computing  a  locally  optimal  quantizer  for  a 
distributed  estimation  system.  Briefly,  one  starts  with  an 
arbitrary  initial  quantizer  and  computes  a  function  /i(1) 
given  by  (3).  Using  and  the  initial  partition,  a  new 
partition  is  generated  using  (4).  One  then  starts  over  and 
uses  the  new  partition  to  generate  according  to  (3), 
and  so  on.  The  algorithm  stops  when  the  mean-square  er¬ 
ror  given  by  (2)  when  h  —  /»(ml  (m  is  the  iteration  number) 
is  small  enough. 

As  the  preceding  discussion  indicates,  the  computational 
size  of  this  problem  grows  exponentially  with  the  number 
of  sensors  n.  Below  we  impose  constraints  on  h  so  that  the 
size  of  the  problem  of  finding  a  locally  optimal  quantizer 
grows  linearly  with  n. 


E[  \X  -  h(Z\, ....  Zn)|2  ];  (2) 

hence,  the  optimal  h  is  the  conditional  expectation, 

*( - 1 - - )  =  E[-Y  |  Z\  =  ;i - ,Z„ 

(f  we  set  =  ij  +  l . in  =  zn  +  1  and  let 

B  —  .4u,  x  x  .4nt^, 
then  this  conditional  expectation  is  given  by 


V)  =  yt . Yn  =  yn  I  dFy,  y„(yi . yn) 

P( Yi  £  ,4 u, . V„  6  .4nl, ) 


(3) 


III.  Constraining  the  Fusion  Center 


Our  approach  [4]  is  to  constrain  the  computational  capa¬ 
bilities  of  the  fusion  center  a  priori  as  follows.  We  require 
that 

fl 

X  =  h(Zx,...,Zn)  =  ]T>(Zt),  (5) 

k  =  t 


where 


N 

9k(Zk)  ±  £cki/{i_1}(Zt). 

i  =  t 


(6) 


Remark:  In  spite  of  the  sums  in  (5)  and  (6),  the  fusion 
center  is  performing  a  nonlinear  operation  on  the  input 


data  Zx . Zn  In  fact,  since  the  Zk  are  discrete  random 

variables,  the  set  of  possible  inputs  does  not  constitute  a 
vector  space  over  IR.  Similarly,  each  g *  in  (6)  is  a  nonlinear 
function  of  Zk,  and  in  (1),  Zk  >8  a  nonlinear  function  of  V*. 


Recalling  that  Z\  is  a  function  of  Y|  (cf.  (1)),  we  can  use 
the  smoothing  property  of  conditional  expectation  to  write 

J,  =  E[g,(Z,){gi(Zl)-2(ElX\Y,]-YlZ[9k(Zk)\Y,])}). 

k*l 


Combining  (5)  and  (6),  and  recalling  that  Zk  —  i  —  1  if 
and  only  if  Y*  G  At,,  we  have 

X  =  (7) 

k-\  i=i 


Clearly,  A'  is  a  nonlinear  function  of  Y\, . . . ,  Yn.  However, 
if  the  partitions  at  the  sensors  are  fixed,  choosing  the  {c**} 
that  minimize  E[|.Y  -  A' | 2 ]  is  a  linear-estimation  problem 
whose  solution  is  given  by  the  usual  normal  equations.  In 
this  case,  we  will  have  ,Vn  equations  in  ;Vn  unknowns. 
Hence,  the  number  of  equations  will  grow  linearly  with 
the  number  of  sensors  n  The  moments  needed  to  write 
down  the  normal  equations  are 


ElCn.O'*)/*, .(Y,)] 


P(  Yt  €  .-U,)6,j , 

P(Vt  6  Aki.Y,  G  A,}), 


whpre  =  1  if  i  =  j  and  6t)  =  0  otherwise,  and 


I  =  k , 
(8) 


=  / 
JAk, 


Yk  =  y]dFYk(y)- 


Note  that  one  needs  only  the  joint  distributions  of  the 
form  Fyty,  and  FXYi  and  not  the  entire  joint  distribution 

Fxy, 

Definition  1:  Given  a  partition  for  each  sensor 

k.  we  write 


A  -  (  Ml,  KV=p  •  -  {Ani}fLi)- 


We  denote  the  procedure  of  solving  the  above-mentioned 
linear  estimation  problem  by  SN(  A).  Letting  C  denote  the 
n  x  ,V  matrix  with  elements  ct,,  we  write  C  =  SN(A). 


We  now  discuss  how  to  find  a  good  partition.  If  the 
matrix  C  is  fixed,  it  is  now  very  natural  to  ask  how  the 
best  partition  is  characterized.  To  obtain  the  answer  to 
this  question,  fix  any  /  =  1 , . . . ,  n,  and  write 


We  can  then  write 

Ji  =  V  /  <pn{y)dFY,(y), 

Ja>‘ 

where  fu(y)  =  c,,(ci,  -  2r,(y))  and 

n(y)  =  E[,Y  |  Vi  =  y]  -  E I  Yi  =  y] 

Ml 

=  E[  A  |  Y,  =  y] 

N 

-ZEc*>p(Y‘eAwiy'®y)-  (11) 

Ml  ;  =  i 

Clearly,  if  the  {c^}  are  fixed  for  all  k  and  i,  and  if  the 
partitions  {,4tl}£L,  are  fixed  for  k  96  /,  then  we  should  put 

y  G  An  <=>  <pu(y)  <  Vli'(y)  for  all  »'  =  1 - N. 

If  we  assume  that  cji  <  ••  <  cis,  then  this  is  equivalent 

to 

a  f  _  C|,_i+C/j  ,  .  Cli  +  C/,i  +  i  \ 

An  =  jyGlR:  -  <  n(y)  <  - - - (12) 

(The  choice  of  <  and  <  is  arbitrary  and  is  made  so  that 
the  will  be  disjoint.)  Observe  that  the  function 

rj  depends  on  the  for  all  k  ^  Also,  the  set 

,4(,  in  (12)  is  not  an  interval,  but  rather  the  inverse  im¬ 
age  of  an  interval.  It  is  also  important  to  observe  that 
to  compute  n  for  /  =  1, . . .  ,n  only  requires  knowing  the 
two-dimensional  joint  distributions  Fxy,  an<I  FYkY,  for  all 
k  and  /.  An  important  consequence  of  this  fact  is  that 
if  we  decide  to  add  another  sensor  to  measure,  say  Yn  +  i, 
our  prior  knowledge  of  Fxy,  and  FYky,  for  k,l  <  n  can  be 
reused.  Of  course,  we  would  still  need  to  obtain  Fxy,+, 
and  FYkYn+,  for  k  =  1, . . . ,  n. 

We  conclude  this  section  with  a  final  definition  and  a 
remark. 


E [ ( A'  -  A' | 2 ]  =  E[|.Y-y,(Z,)-£y*(Zt)|2].  (10) 

Ml 

The  right-hand  side  of  (10)  can  be  expanded  into  nine 
terms:  however,  only  five  terms  will  involve  Zi .  Denoting 
the  sum  of  these  five  terms  by  Ji,  we  have 

J,  =  E[ffi(Z,){jr,(Z,)-2(.Y  -]r<MZ*))}l. 

Ml 


Definition  2:  We  introduce  the  procedure  Ui(C,  A).  Re¬ 
call  our  notation  in  Definition  1.  Let  A  =  Ui(C,  A)  be  ob¬ 
tained  from  A  by  replacing  {A/jJfLj  with  ,  where 

each  An  is  given  by  the  right-hand  side  of  (12). 

Remark:  If  n  =  1  and  AT  =  Yi,  then  the  normal  equa¬ 
tions  reduce  to 

P  (Yx€Au)cx,  =  E[Y,h,.(Yi)]. 


5" 


In  other  words,  we  recover  the  classical  Lloyd-Max  condi¬ 
tions  for  locally  optimal  quantizers  [6,  7], 


where  X,  W\,  and  W?  are  statistically  independent  were 
considered. 

Example  (8,  Example  8]:  Let  X  have  density 


■  I ? 


-  cos(^)], 


kl  < 

otherwise, 


where  d  ss  0.3419  is  a  normalization  constant  and  6  =  2 
(see  Fig.  2).  We  let  Wi  and  W 2  have  the  same  density 
except  that  6=1. 


IV.  The  Design  Algorithm 

Using  the  basic  procedures  SN  and  U\,  l  —  1, . . . ,  n,  de¬ 
fined  in  the  preceding  section,  there  are  two,  almost  iden¬ 
tical,  algorithms  for  generating  approximately  locally  op¬ 
timal  quantizers  for  distributed  estimation  systems. 

Algorithm  1 

Let  A  =  ({AuJfLj, . .  ,  {.4„,},^1)  be  given. 
C  :=  SN(A) 
loopl:  for  /  =  1  to  n 

A  :=  U,(C,  A) 
next  l 
C  :=  SN(A) 

if  stopping  criterion  not  met,  go  to  loopl 
end 


Algorithm  2 

Let  A  =  ({Aw}",, . . . ,  be  given. 

C  :  =  SN(A) 
loop2:  for  /  =  1  to  n 

A  :=  Ui(C,  A) 

C  :  =  SN(A) 

next  / 

if  stopping  criterion  not  met,  go  to  loop2 

end 

While  the  preceding  algorithms  appear  simple  enough, 
their  implementation  is  nontrivial.  The  two  main  difficul¬ 
ties  in  implementing  the  algorithms  are  the  computation 
of  the  function  ri(y)  in  (11)  and  the  characterization  of  the 
inverse  images  in  (12).  Note  that  even  if  X,  Y\, . . . ,  Yn  are 
jointly  Gaussian,  we  cannot  write  (11)  in  closed  form  even 
if  the  Akj  are  intervals.  Hence,  the  sets  A»,  in  (12)  must 
be  determined  numerically  and  then  a  description  of  them 
must  be  stored  in  a  suitable  data  structure. 

A  set  of  programs  has  been  developed  [8]  to  implement 
Algorithms  1  and  2  when  provided  with  subroutines  to 
compute  the  particular  moments  and  probabilities  for  a 
given  situation.  Several  examples  of  the  form 

n  =  A  +  Wk,  k  =  1,2, 


With  N  =  8  (3-bit  quantizers),  the  Lloyd-Max  partitions 
for  the  random  variables  Y*,  k  =  1,2,  are  identical  and  are 
given  by 

Atl  =  (-00,-1.3273] 

At2  =  (-1.3273,-1.2419] 

Ai3  =  (-1-2419,-1.0374] 

Ai4  =  (-1.0374,0] 

A*s  =  (0,1.0374]  ’ 

Afc6  =  (1.0374,1.2419] 

At7  =  (1.2419,1.3273] 

Am  =  (1.3273,oo). 

Solving  the  normal  equations  for  C  —  SN(A),  we  have 
E[|A-X|J]  =  0.18129. 

Note  that  the  minimum  mean  square  error  achievable  by 
a  linear  estimator  is  0.18534.  Thus,  just  using  the  Lloyd- 
Max  quantizers  and  doing  linear  estimation  on  Zi,Zi  can 
do  better  than  pure  linear  estimation.  Using  Algorithm  1, 
we  initialized  A  to  the  Lloyd-Max  partition  in  (14).  After 
5  passes  through  the  loop  in  Algorithm  1,  the  partitions 


I 


were 

An  =  (-00,  -2.1802] 

.4 j 2  =  (-2.8102, -1.7950]U  (-1.3697, -0.8177] 
4i3  =  (-1.7950, -1.3697]U  (—0.8177,  -0.4109] 
.4U  =  (-0.4109, 0.0005686] 

415  =  (0.0005686,0.4130] 

.4 1 6  =  (0.4130, 0.8200]  U  (1.3739, 1.7979] 

.4  it  =  (0.8200, 1.3739]  U(  1.7979, 2. 1834] 

.4  is  =  (2.1834,oo) 

and 


.421 

= 

(-oo.-2.1056] 

-422 

= 

(  —  2. 1056,  —0.6504] 

■423 

= 

(-0.6504,  -0.3207] 

-424 

= 

( —0.3207.  —0. 1805] 

.425 

= 

(—0  1805.  0  2863] 

.426 

= 

(0.2863,0.0240] 

.427 

= 

(0.6240,2.0960] 

•426 

— 

(2.0960,  oo). 

The  minimum  mean  square  error  for  these  partitions  is 
E[  |.Y  -  A’j2 ]  =  0.12655, 

which  is  more  than  a  30%  improvement  over  the  perfor¬ 
mance  of  the  Lloyd-Max  partition  and  over  pure  linear 
estimation. 

Remarks:  (i)  After  5  passes  through  the  algorithm,  the 
mean  square  error  was  not  significantly  reduced. 

(ii)  The  final  partitions  for  the  sensors  are  not  the  same, 
even  though  sensors  1  and  2  play  interchangeable  roles  in 
this  example.  The  reason  for  this  is  that  the  algorithm 
treats  one  sensor  at  a  time. 

(i»«)  As  a  general  rule,  it  was  found  in  [8]  that  Algo¬ 
rithm  2  yielded  results  almost  identical  to  those  of  Algo¬ 
rithm  1. 
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V.  Conclusion 

We  have  developed  an  algorithm  for  the  design  of  a  dis¬ 
tributed  estimation  system  with  n  sensors  and  a  single  fu¬ 
sion  center  that  is  subject  to  communication  and  computa¬ 
tion  constraints.  The  algorithm  uses  only  bivariate  prob¬ 
ability  distributions  and  yields  locally  optimal  estimators 
that  satisfy  the  required  system  constraints. 

While  this  work  was  motivated  by  problems  in  sensor 
fusion,  the  ideas  can  also  be  applied  in  a  general  nonlin¬ 
ear  estimation  context.  In  other  words,  estimators  of  the 
form  (7)  constitute  a  class  of  nonlinear  estimators,  and  the 
algorithm  presented  here  can  be  used  to  obtain  a  locally 
optimal  nonlinear  estimator  from  this  class. 
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I 


1.  Introduction 


Distributed  estimation  systems  lie  at  the  heart  of  many  important  applications.  For 
example,  distributed  estimation  systems  are  employed  by  public  power  utilities  to  determine 
various  loads  and  voltages  both  inside  and  outside  their  services  areas.  Improved  estimation 
systems  would  permit  them  to  provide  electricity  more  economically.  Another  application  of 
distributed  systems  is  in  the  area  of  air  travel  control.  In  this  case,  many  dispersed  radars 
monitor  aircraft  in  order  to  provide  guidance  and  to  avoid  collisions.  Better  estimation 
systems  would  provide  safer  air  travel. 

In  this  report,  we  consider  a  distributed  estimation  system  consisting  of  n  sensors  that 
transmit  their  findings  to  a  common  fusion  center  over  a  channel  of  positive,  finite  capacity. 
Using  such  a  system,  it  is  desired  to  estimate  some  quantity  X,  which  cannot  itself  be 
observed.  Because  the  communication  channel  has  finite  capacity,  it  is  necessary  for  the 
sensors  to  quantize  their  measurements  before  transmission. 

After  some  preliminary  definitions  in  Section  2,  we  give  a  detailed  overview  of  the  con¬ 
strained  estimator  approach  in  Section  3.  In  Section  4  we  present  an  algorithm  for  quan¬ 
tization  in  distributed  estimation  systems  (the  multi-sensor  case)  for  which  the  theory  was 
outlined  in  Section  3.  Section  4  also  contains  some  numerical  results  obtained  using  this 
algorithm.  Finally  in  Section  5  we  give  an  overview  of  the  code  used  to  implement  the 
algorithm  presented  in  Section  4. 

2.  Background  and  Notation 

Consider  the  distributed  estimation  system  shown  in  Figure  1.  This  system  consists  of 
n  sensor  platforms  that  transmit  their  findings  over  a  communication  channel  to  a  common 
fusion  center.  The  sensors  do  not  communicate  with  each  other.  Using  this  system  some 
quant  ity,  X ,  is  to  be  estimated.  However,  only  the  sensor  measurements  Y\, . . . ,  Yn,  which  are 
related  to  X ,  are  available.  Suppose  that  the  findings,  Y\, . . . ,  Yn,  of  the  n  sensor  platforms 
are  transmitted  to  the  fusion  center  over  channels  capable  of  transmitting  at  most  R  bits  per 
second  during  a  finite  time  interval  [0,  T].  Then  each  sensor  k  communicates  the  real  number 
Yk  by  sending  RT- bit  word  to  the  fusion  center.  Hence,  it  is  necessary  for  the  sensors  to 


quantize  their  measurements. 


Y\  K 


Estimate  of  X 


Figure  1:  A  Distributed  Estimation  System. 

Let  X,  Y\ , . . . ,  Yn  be  real- valued  random  variables  on  some  probability  space  (Q,P,  P). 
Since  each  measurement  V*  must  be  mapped  into  an  RT- bit  word,  it  is  convenient  to  set 
N  =  2RT .  We  assume  that  each  sensor  k  is  equipped  with  a  partition,  {A*, •}£•!,  and  that 
sensor  k  transmits  the  random  variable 

2*  =  Bi-lJWtt)  (1) 

1=1 

to  the  fusion  center  without  error.  (Here  Ia{v)  denotes  the  indicator  function  of  the  set  A ; 
i.e.,  />t(y)  =  1  if  y  €  A  and  Ia(v)  =  0  otherwise.)  Note  that  the  event  {Zk  =  *  —  1}  is  the 
same  as  the  event  {H  €  A*,}. 

Below  we  give  an  overview  of  the  constrained  estimator  approach,  which  yields  an  algo¬ 
rithm  for  quantization  in  distributed  estimation  systems.  We  should  note  that  the  goad  of 
quantization  in  such  systems  is  to  provide  the  best  estimate  of  the  unobservable  X  given 
the  measurements  Yu . . . ,  Yn  related  to  X.  Our  goal  is  not  to  reconstruct  the  measurements 
Yi, . . . ,  Yn  at  the  fusion  center. 


n 


3.  The  Constrained  Estimator 


Our  approach  is  to  constrain  the  computational  capabilities  of  the  fusion  center  a  priori 
as  follows.  Denote  the  fusion  center  output  by  X.  We  require  that 


x  =  £s*(z»). 

i=l 


where 


N 


9k(Zk )  =  £cfct/{,_1}(Z*). 


i=i 


In  other  words,  since  Zk  =  i  —  1  if  and  only  if  Yk  £  Aki, 

x  =  EEcaun). 

/fc=l  »  =  1 


(2) 


Clearly  X  is  a  nonlinear  function  of  Yk.  It  is  convenient  to  rewrite  (2)  in  matrix  notation. 
Let  the  matrices  h  and  a  be  defined  as: 

h  &  i iA„m,ujY, . n.„(K)]T 

a,  —  [  Cji  ,  C\2  ,  •  •  •  ,  Cnl  i  •  •  •  j  &nN  ]  • 


Then  the  fusion  center  output  becomes, 


X  =  aTh. 


(3) 


When  the  partitions  at  the  sensors  are  fixed,  the  problem  of  choosing  a  to  minimize 

E[  \X  —  X\2  ]  is  a  simple  linear  estimation  problem  whose  solution  is  given  by  the  normal 

equations.  That  is,  the  optimal  a  is  the  solution  of 


Ra  =  r, 


(4) 


where  R  =  E[hhT]  and  r  =  E[A"ft].  To  solve  the  normal  equations,  we  have  Nn  equations 
in  Nn  unknowns.  Hence,  the  number  of  equations  will  grow  linearly  with  the  number  of 
sensors  n.  For  example  if  n  =  4  and  N  =  2,  there  will  be  8  equations.  The  moments  needed 
to  write  down  the  normal  equations  are 


’  P(Yk  €  Akt)6t],  \U  =  k, 

P(Yk  €  Aki,Yt  £  Au),  if  l^k, 


(5) 


where  Sij  =  1  if  i  =  j  and  0  otherwise,  and 

E[X/,..(tt)I  =  /  t[X\Yt  =  y)dFr.(y).  (6) 

JAk, 

Definition  1.  Given  a  partition  for  each  sensor  k ,  we  write 

A  =  ({An}?,,. 

We  denote  the  procedure  of  solving  the  above-mentioned  linear  estimation  problem  by 
SN(A).  Letting  C  denote  the  n  x  N  matrix  with  elements  cjt,,  we  write  C  =  SN(A). 
Clearly  C  and  a  contain  the  same  data,  simply  rearranged. 

How  to  find  a  good  partition.  If  the  matrix  C  is  fixed,  it  is  now  very  natural  to  ask  how  the 
best  partition  is  characterized.  To  obtain  the  answer  to  this  question,  fix  any  l  =  I, . . .  ,n, 
and  write 

E(|.V-A-|!]  =  E[|.Y-!,,(ZI)-5>(Z0Ij].  (7) 

k*l 

The  right-hand  side  of  (7)  can  be  expanded  into  nine  terms;  however,  only  five  terms  will 
involve  Z{.  Denoting  the  sum  of  those  five  terms  by  J/,  we  have 

Ji  =  E[J,(Z,){„(Z|)-2(X-  5>(Z*))}]. 

MI 

Recalling  that  Zi  is  a  function  of  yj  (cf.  (1)),  we  can  use  the  smoothing  property  of  conditional 
expectation  to  write, 


J,  =  E[*(Z,){J,(Z,)-2(E[X|y,]-j;E[w(Z,)l«])}]- 

Ml 

We  can  then  write 

Ji  =  £/,  <*{.*-2(E[X|«  =  »|-EE[s*<Z‘)IH  =  vl)}  dFr,(y). 

1=1  JAl'  k*l 

' - - - - 

Clearly,  if  the  {c^, }  are  fixed  for  all  k  and  i,  and  if  the  partitions  are  fixed  for  k  ^  /, 

then  we  should  put 


y  €  An 


<Pi i(y)  <  <Pii>(y)  for  all  t'  = 


(8) 


We  can  obtain  a  simpler  expression  for  (8)  by  setting 


n(y)  =  £[X\Yl  =  y]-YdZ.[9k{Zk)\Yl  =  y] 

k*l 

=  E[X\Yl  =  y]-'£'tckjP(YkeAkj\Yl  =  y). 

k*l  ;  =  1 


(9) 


We  can  then  write  =  c/,(c/,  —  2r/(y)).  If  we  assume  that  cn  <  ■  ■  •  <  cim,  then  (8)  is 

equivalent  to 


An  =  jy  6  IR  : 


Cl.i- 1  +  Cli  Cli  +  C/,,+ 1 

■  <  n(y)  <  ■ 


2  —  2  y  (10) 

Observe  that  the  function  r/  depends  on  the  {cjtj}^.,  for  all  k  l .  Also,  the  set  Ai,  in  (10) 

is  not  an  interval,  but  rather  the  inverse  image  of  an  interval. 


Definition  2.  Recall  our  notation  in  Definition  1.  Let 

me,  a)  *  <{*„>£, . <^)£..  (-W&, . <*-}£,), 

where  each  An  is  given  by  the  right-hand  side  of  (10). 

Finally  to  evaluate  the  performance  of  the  constrained  estimator  algorithm  the  mean 
square  error  is  computed.  If  X  satisfies  (3)  and  a  satisfies  (4),  the  mean  square  error  is: 

E[|X-X|2]  =  E[X2  +  X2  —  2XX] 

=  E[X2  +  (aTh)2  -2aTXh] 

=  E[  X2  ]  —  2arr  +  aTRa 


=  E[  X2  ]  —  arr 

=  E[xji-EEcitE[x/4(n)|. 

fc=i  >=i 

Remark.  If  n  =  1  and  X  =  Y\,  then  the  normal  equations  reduce  to 

P(Y  €Au)cu  =  E WAuiYi)], 


(ID 


or 


ci,  =  ^ 


[  ydFvAy) 

J  At. 


P(*i  €  Au) 


w 


Further,  ri(y)  =  y,  and  so 


A  f  /-  TD  +  C>*  ^  ^  CU  +  Cl,i+l  | 

Au  =  |S/  €  1R  :  - - -  <y<  - - - j. 

In  other  words,  we  recover  the  classical  Lloyd-Max  conditions  for  locally  optimal  quantizers 

[1,2]. 

4.  Numerical  Results 

Using  the  basic  procedures  SN  and  Ui,  l  =  l,...,n,  defined  in  the  preceding  section, 
there  is  an  obvious  algorithm  for  obtaining  (approximate)  locally  optimal  quantizers  for 
distributed  estimation  systems. 

Algorithm  1 

Let  A  =  ({i411-}jl1,...,{i4ni}^1)  be  given. 

C  :=  SN(A) 
loopl:  for  /  =  1  to  n 

A  :=  Ui(C,  A) 
next  / 

C  :=  SN(A) 

if  stopping  criterion  not  met,  go  to  loopl 
end 

While  the  preceding  algorithm  appears  simple  enough,  its  implementation  is  not  trivial. 
The  two  main  difficulties  in  implementing  the  algorithm  are  the  computation  of  the  func¬ 
tion  r/(y)  in  (9)  and  the  characterization  of  the  inverse  images  in  (10).  Note  that  even  if 
X,  Yi, . . . ,  Yn  are  jointly  Gaussian,  we  cannot  write  (9)  in  closed  form  even  if  the  Akj  are 
intervals.  Hence,  the  sets  An  in  (10)  must  be  determined  numerically  and  then  a  description 
of  them  must  be  stored  in  some  data  structure. 

Algorithm  1  has  been  programmed,  and  we  have  run  several  numerical  examples.  In  all 
the  examples  below,  we  took 

=  X  +  Wk,  *=l,...,n,  (12) 

where  Wn  are  statistically  independent. 


Example  1.  We  first  took  n  =  1  sensor,  N  —  i  (2-bit  quantizers),  and  Yx  =  X  ~  mi); 
i.e,  X  is  normal  with  mean  0  and  variance  1.  In  other  words,  we  used  our  algorithm  to  solve 
the  classical  Lloyd-Max  problem.  From  several  initial  partitions,  our  algorithm  quickly 
converged  to  the  optimal  quantizer, 

Ak !  =  (-oo,-0.9816] 

Aki  =  (-0.9816,0] 

=  (0.0.9816] 

.4,4  =  (0.9816,oo), 

which  was  originally  reported  in  [2,  Table  I]. 

Example  2.  Suppose  ti  =  2  sensors,  N  =  4  (2-bit  quantizers),  where  X  ~  ./V(0, 0.9)  and 
IVi,  W2  ~  .Y(0, 0.1).  Since  Y\  and  K2  are  iV( 0, 1)  we  take  the  initial  partition  for  each  sensor 
to  be  the  Lloyd-Max  partition  found  in  Example  1.  Taking  C  =  SN(A),  we  find  that 

E[|X-Xp]  =  0.113. 

After  6  passes  through  the  loop  in  Algorithm  1,  we  find  that 

.4,,  =  (-oo,-1.28]  An  =  (-oo,-0.843] 

A 12  =  (-1.28,0]  A22  =  (-0.843,0] 

and 

.4,3  =  (0,1.28]  .4 23  =  (0,0.843] 

.4,4  =  (1.28,cx))  .4  24  =  (0.843,  oo) 

and 

E[\X-X\2}  =  0.107, 
which  is  about  a  5%  reduction  in  the  mean  square  error. 

Example  3.  Let  n  =  2  sensors,  N  =  4  (2-bit  quantizers),  where  this  timeX  ~  U(— 0.9, 0.9); 
i.e.,  uniform  over  the  interval  (—0.9, 0.9).  We  take  Wx,Wi  ~  {/(— 0.11,0.11).  Clearly,  F,, 
k  =  1,2,  has  a  density  given  by  the  convolution  of  two  uniform  densities.  We  initialize  the 


partition  of  each  sensor  to  its  Lloyd-Max  partition,  which  in  this  case  is  ( k  =  1,2) 

Ak\  =  (-oo,-0.4546] 

Ak2  =  (-0.4546,0] 

Aa  =  (0,0.4546] 

Ak\  =  (0.4546,  oo). 

Taking  C  =  SN(A),  we  find  that 

E[\X-X\2]  =  0.01358. 


After  21  passes  through  the  loop  of  Algorithm  1,  we  find  that 


and 


n  =  (-oo,-0.6390] 

A2i  =  (-oo,-0.3929] 

i2  =  (-0.6390,-0.1062] 

and 

A22  =  (-0.3929,-0.1450] 

13  =  (-0.1062.0.6224] 

A23  =  (-0.1450,0.3620] 

i4  =  (0.6224,oo) 

A2  4  =  (0.3620,  oo) 

E[  \X  -  X]2]  =  0.00895, 


which  is  a  reduction  of  about  34.5%  in  the  mean  square  error.  We  should  note  that  after 
only  8  passes  the  mean  square  error  had  been  reduced  by  about  33%. 


Example  4.  We  modify  the  previous  example  by  increasing  the  partition  size  from  N  =  4 
(2-bit  quantizers)  to  iV  =  64  (8-bit  quantizers).  We  initialize  the  partition  of  each  sensor  to 
its  Lloyd-Max  partition.  Taking  C  =  SN(A),  we  find  that  the  initial  mean  square  error  is, 

E[  |A"  —  AT|2]  =  0.00197. 

Algorithm  1  did  not  improve  substantially  the  initial  mean  square  error;  however,  the  value 
0.00197  is  smaller  than  the  minimum  mean  square  error  obtainable  by  linear  estimation, 
which  is  0.0020012. 


Example  5.  Suppose  n  =  2  sensors,  N  =  4  (2-bit  quantizers),  where  X 
Wi  ~  C7( — 0.11,0.11)  and  the  density  for  W2  is 

[  2tan-1(6a)(l  -I-  {aw)2)  ’  ^1  —  ^’ 


fwAw)  =  s 


~  £/(— 0.9, 0.9), 


(13) 


0, 


otherwise, 


where  b  =  0.11  and  a  =  1.  We  initialize  the  partition  of  each  sensor  to  its  Lloyd-Max 
partition.  We  have 


An  =  (-00,-0.4546] 

A2  i  — 

(—oo,  -0.4456] 

An  =  (-0.4546,0] 

A  22  = 

(-0.4456,0] 

and 

A, 3  =  (0,0.4546] 

A  23  = 

(0,0.4456] 

Ah  =  (0.4546,  oo) 

A24  = 

(0.4456,  oo). 

Taking  C  =  SN(A),  we  find  that 

E[|* 

-X\2}  =  0.0136. 

After  21  passes  through  the  loop  in  Algorithm  1,  we  find  that 

An  =  (—oo,—0.6390] 

A2l  — 

(-oo,-0.3927] 

Au  =  (-0.6390,0.1066] 

a22  = 

and 

(-0.3927,-0.1447] 

An  =  (0.1062,0.6226] 

A23  = 

(-0.1447,0.3624] 

Au  =  (0.6226,  oo) 

K> 

A. 

II 

(0.3624, oo) 

and 

E[ \X  ■ 

-X\2}  =  0.00894, 

which  is  about  a  34.3%  reduction  in  the  mean  square  error.  We  should  note  that  after  only 
8  passes  the  mean  square  error  was  reduced  by  33.5%.  These  results  are  almost  identical 
to  those  in  Example  3;  this  is  not  surprising,  since  taking  a  =  1  in  (13)  implies  that  fw3  is 
nearly  flat,  and  hence  W?  is  nearly  uniform  as  in  Example  3. 

Example  6.  Same  as  Example  5,  except  we  replace  a  =  1  by  a  =  32.  We  initialize  the 
partition  of  each  sensor  to  the  Lloyd-Max  partition,  which  is 


An  =  (—oo,—0.4546] 

A2i  =  (-oo,-0.4481] 

A,  2  =  (-0.4546,0] 

and 

A22  =  (-0.4481,0] 

A, 3  =  (0,0.4546] 

A 23  =  (0,0.4481] 

Ah  =  (0.4546,  oo) 

A24  —  (0.4481,  oo). 

Taking  C  =  SN(A),  we  find  that 


E[\X-X\2}  =  0.0136. 


After  10  passes  through  the  loop  in  Algorithm  1,  we  find  that 


An  =  (-oo,-0.6374] 

A2\  =  (— oo,  —0.3852] 

A, 2  =  (-0.6374,0.1462] 

and 

A22  =  (-0.3852,0.0913] 

Al3  =  (0.1462,0.6160] 

A23  =  (0.0913,0.3496] 

Ah  =  (0.6160,  oo) 

A2 4  =  (0.3496,  oo), 

and 

E[|X-X|2]  =  0.00797, 


which  is  about  a  41.4%  reduction  in  the  mean  square  error.  We  should  note  that  after  only 
2  passes  the  mean  square  error  was  reduced  by  37.2%. 


After  15  passes  through  the  loop  in  Algorithm  1,  we  find  that 


An 
A\ 2 

Al3 
A 14 


(-oo,-0.6156] 

(-0.6156,-0.1088] 

and 

(-0.1088,0.6234] 

(0.6234,  oo) 


A2 !  =  (-00,-0.3570] 
A22  =  (-0.3570,0.1328] 
A23  =  (0.1328,0.3721] 
A24  =  (0.3721,  oo), 


■  o 

■  / 


and 


E[|A-A|2]  =  0.00938, 

which  is  about  a  31%  reduction  in  the  mean  square  error.  We  should  note  that  after  3  passes 
the  mean  square  error  was  reduced  by  25.5%. 


Example  8.  Let  n  =  2  sensors,  N  =  8  (3-bit  quantizers),  where  this  time  A"  has  the  density 
function  shown  in  (14)  where  b  =  2.  We  take  W,,  W2  to  have  the  same  density  function  of  A, 
except  we  replace  b  =  2  by  6  =  1.  We  initialize  the  partition  of  each  sensor  to  its  Lloyd-Max 
partition,  which  in  this  case  is  (k  =  1,2) 

.4a  =  (-oo.-l.3273] 

.4*2  =  (-1.3273,-1.2419] 

.4«  =  (-1.2419,-1.0374] 

.4*4  =  (-1.0374,0] 

.4*5  =  (0,1.0374] 

.4*6  =  (1.0374,1.2419] 

.4*7  =  (1.2419,1.3273] 

/Us  =  (1.3273,oo). 

Taking  C  =  SN(A),  we  find  that 

E[ | A  -  A]2]  =  0.18129. 


After  only  5  passes  through  the  loop  of  Algorithm  1,  we  find  that 


An  =  (-oo.-2.1802] 

Au  =  (-2.8102,-1.7950]  U  (-1.3697,-0.8177] 
Ai3  =  (-1.7950, -1.3697]  U  (-0.8177, -0.4109] 
Au  =  (-0.4109.0.0005686] 

.4,5  =  (0.0005686,0.4130] 

A, 6  =  (0.4130,0.8200]  U  (1.3739,1.7979] 

A, 7  =  (0.8200, 1.3739]  U  (1.7979,2.1834] 

A,s  =  (2.1834,oo) 


A21  = 

(-oo,-2.1056] 

A22  = 

(-2.1056,-0.6504) 

•4  23  = 

(-0.6504,-0.3207) 

A24  = 

(-0.3207,-0.1805) 

A25  = 

(-0.1805,0.2863) 

•4  26  = 

(0.2863,0.6240) 

■427  = 

(0.6240,2.0960) 

A28  = 

(2.0960, 00) 

E[|-V  - 

-X\2}  =  0.12655, 

which  is  a  reduction  of  about  30.1%  in  the  mean  square  error.  We  should  also  point  that  the 
value  0.12655  is  a  31.7%  reduction  of  the  minimum  mean  square  error  obtainable  by  linear 
estimation,  which  is  0.18534. 

5.  Conclusion 

The  best  estimate  of  some  quantity  X,  in  a  distributed  estimation  system,  given  some 
measurements  say  Y\, ...  ,Yn,  is  the  conditional  expectation  of  X  given  Y\, . . . ,  Yn,  denoted 

E[  A’IVi, - Yn].  However,  the  functional  form  of  E[  X\Y\, . . . ,  Yn  ]  as  a  function  of  Y\, . . .  ,Yn 

is  often  difficult,  if  not  impossible  to  determine  and  it  is  not  possible  to  transmit  real- valued 
quantities  without  distortion.  In  this  report  we  solved  those  problems  by  showing  the  feasi¬ 
bility  and  efficiency  of  an  algorithm  that  was  specifically  directed  to  address  the  distributed 
nature  of  the  system.  In  Section  3  we  presented  the  constrained  estimator  approach,  which 
is  an  algorithm  for  quantization  for  distributed  estimation  systems.  The  numerical  results 
in  Section  4  indicated  that  in  general  the  constrained  estimator  performs  better  than  the 
Lloyd-Max  estimator  when  N  is  small  (N  =  4,  2-bit  quantizers)  as  shown  in  Examples  2,3, 
and  5-7,  however  when  N  is  large  (N  =  64,  8-bit  quantizers)  the  constrained  estimator  does 
not  improve  the  results  of  the  Lloyd-Max  estimator  significantly  as  shown  in  Example  4. 

An  other  conclusion  can  be  drawn  from  the  numerical  results.  The  constrained  estimator 
is  a  nonlinear  estimator  as  pointed  out  in  Section  2.  This  is  shown  in  Example  4  and  more 


clearly  in  Example  8.  In  fact,  the  nonlinear  constrained  estimator  performs  better  or  in  the 
worst  case  as  well  as  the  linear  estimation  provided  that  N  is  chosen  appropriately. 


Appendix  A.  An  Overview  of  the  Code 

In  this  section  we  discuss  in  detail  the  code  that  was  used  to  implement  the  constrained 
estimator  algorithm  presented  in  Section  4.  But  first  we  turn  our  attention  to  the  few 
changes  that  we  made  to  the  equations  presented  in  Sections  3  and  4.  In  view  of  the 
structure  described  in  (12),  it  is  convenient  to  express  the  quantities  in  (5),  (6),  and  (9)  in 
terms  of  the  univariate  densities  fx(x)  and  fw,{w).  We  find  that 

P{Yk  €  A)  =  f  P{Wk  €  A  -  x)fx(x)dx,  (15) 

J -OO 


P(V;  €  A,Y,  €  B)  =  P(Wke  A-x)P(W,e  B-x)fx{x)dx,  k  ±  l, 

J  —  OO 

/OO 

P(Wk€A-x)fx(x)dx, 

-OO 


/  P(Wfc  €  A  -  x)fWl(y  -  x)fx{x)dx 
P(n  e  A\Y,  =  y)  =  ^=22 - - - , 


fr\y) 


E[X\Yl  =  y)  = 


roo 

/  x  fw,{y  -  x)  fx{x)dx 

J -OO 


fy,(y) 


where,  of  course, 


and 


roo 

JY,{y)  =  /  fw,(y  ~  x)  fx{x)dx, 

J — OO 


P{WkeA-x )  =  [  fwk(w)  dw. 

J  A—r 

In  particular,  if  A  =  (a,  6],  A  —  x  =  (a  —  x,  6  —  x],  and 

rb-x 

P(Wk  €A-x)  =  /  fwk{w)dw. 

J a-r 


(16) 

(17) 

(18) 

(19) 

(20) 

(21) 


Now  that  we  mentioned  the  few  changes  needed  to  make  the  implementation  of  our 
algorithm  easier,  we  are  ready  to  discuss  the  code  used  in  detail.  As  shown  in  Algorithm  1, 


the  first  element  that  is  needed  to  implement  this  algorithm  is  an  arbitrary  initial  partition 
A.  The  partition  A  is  defined  as: 

A  =  ({*„}£, . {a„,)£,) 

where  {A*,}^x  is  a  partition  for  sensor  k.  We  call  Aki  the  ith  subpartition  of  sensor  k. 

We  now  discuss  the  procedure  that  read  and  stored  the  initial  partition  A.  The  initial 
partition  is  stored  in  the  file  start.  The  lower  and  upper  limits  of  the  intervals  of  each 
subpartition  1  through  N  are  listed  in  succession  and  are  separated  by  the  pair  -1,-1.  This 
is  repeated  for  each  sensor  1  through  n.  We  should  note  that  the  interval  (—00,00)  is  stored 
as  the  pair  0,-1  while  the  interval  (  —  00,  a)  is  stored  as  the  pair  a +  2,  a,  and  (a,  00)  is  stored 
as  the  pair  a,  a  —  3. 

For  example,  if  n  =  1  and  N  =  2  then  the  partition  A  is  defined  as  A  =  ({AnJ^Lj), 
and  if  An  =  (—00, 3]  U  (4,  00)  and  AX2  =  (3,4],  then  the  file  start  will  contain  the 
following  lines: 

5,  3 

4,  1 

-1,  -1 
3,  4 

-1,  "I 

This  listing,  available  in  file  start  is  read  and  then  stored  in  a  3-dimensional  array  AR(I , J ,K) 
by  the  routine  ReadAR.  The  first  index  I  of  the  array  AR(I,J,K)  indicates  which  of  the 
sensors  from  1  through  n  we  are  dealing  with,  the  second  index  J  (J=l,. . .  ,N)  indicates 
which  subpartition  we  are  pointing  to,  and  finally,  the  third  index  K  indicates  if  the  number 
stored  is  the  upper  or  lower  endpoint  of  an  interval  of  a  subpartition  at  J  of  sensor  I.  For 
example,  the  set  An  described  above  in  the  file  start  will  be  stored  in  AR  as  follows: 

AR(  1,1,1)  =  5  AR(1, 1,3)  =  4 

and 

AR(1, 1,2)  =  3  AR(1, 1,4)  =  1, 

and  An  will  be  stored  in  AR  as  follows: 

AR(  1,2,1)  =  3 

AR(1,2,2)  =  4. 


The  routine  ReadAR  reads  the  file  start  and  returns  the  array  AR  containing  the  initial 
partition.  The  key  for  the  proper  operation  of  the  routine  is  the  knowledge  of  the  number 
of  sensors,  n,  denoted  by  NS  in  the  programs,  and  knowledge  of  the  number  of  partitions, 
N  denoted  by  N  in  the  programs.  The  routine  reads  every  pair  of  elements  from  start  and 
stores  them  in  AR(I,J,K)  and  AR(I,J,K+1),  this  process  is  repeated  until  a  pair  of  —Is 
occurs  indicating  the  subpartition  is  updated  and  J  becomes  J+l  and  the  next  subpartition 
is  read.  When  all  subpartitions  of  intervals  of  a  sensor  I  are  read  ( J=N)  and  stored  in  AR,  the 
index  indicating  the  number  of  sensors  is  updated,  and  I  becomes  1+1.  The  whole  process 
is  then  repeated  until  I=NS;  at  this  point  the  partition  A  has  been  stored  in  AR. 

Now  that  the  initial  partition  A  is  read  and  stored  in  the  matrix  AR(I,J,K),  the  next 
step  is  to  compute  the  initial  elements  {cit,}.  To  do  that  we  generate  the  two  matrices  R 
and  r  as  defined  in  (4).  In  the  programs,  R  is  called  R  and  r  is  called  Mean. 

The  routine  GenR(R)  returns  the  matrix  R  defined  in  (4)  using  the  routine  SumP(L,  J  ,K, I) 
that  generates  the  moments  shown  in  (5),  (15),  and  (16).  Two  nested  do  loops  were  used; 
the  first  loop  generated  all  the  elements  of  R  and  stored  them  in  the  N2n 2  x  1  matrix  Temp2, 
the  second  do  loop  stored  appropriately  the  elements  of  Temp2  in  the  Nn  x  Nn  matrix  R. 

The  routine  GenM(Mean)  returns  the  matrix  r  denoted  as  Mean  using  the  routine  SumM 
(K,I)  that  generates  the  moments  shown  in  (17).  In  this  case  one  nested  loop  is  used,  this 
nested  do  loop  will  generate  the  iVn  x  1  matrix  Mean. 

Now  that  the  matrices  R  and  Mean  are  available,  the  routine  Solve(R,Mean,C,MSE) 
computes  the  coefficients  { c*. }  stored  in  the  n  x  N  matrix  C,  where  each  row  indicates 
a  different  sensor  and  the  elements  of  a  given  row  k  are  the  coefficients  of  the  sensor  k 
that  minimizes  E[  \X  —  X|2].  This  routine  also  returns  the  mean  square  error  MSE,  which 
is  computed  using  (11).  Since  the  normal  equations  used  for  obtaining  the  {c*;}  are  ill 
conditioned,  we  employed  the  singular-value  decomposition  (SVD)  as  discussed  in  [3].  To 
obtain  the  SVD  of  R  we  used  the  NAG  routine  F02WEF. 

Finally,  the  routine  SN(AR,C)  will  return  the  n  x  N  matrix  C  containing  the  different 
elements  {c*,}  obtained  by  solving  the  mentioned  linear  estimation  problem  defined  in  (4) 
given  the  initial  partition  A  stored  in  the  array  AR.  This  routine  is  used  in  the  main  program 
and  is  equivalent  to  the  procedure  SN  in  Definition  1. 


We  have  just  shown  how  to  find  the  elements  {c^}  that  will  minimize  E[  \X  —  X\ 2  ]  given 
a  partition  A  using  the  procedure  SN(AR,C).  The  next  step  is  to  present  the  procedures  that 
we  used  to  improve  a  partition  A.  As  we  seek  to  improve  a  partition  A  given  a  matrix  C, 
we  have  a  major  hurdle  to  overcome.  The  first  part  of  that  hurdle  is  the  computation  of  the 
function  n(?/)  in  (9),  the  second  part  is  to  characterize  the  inverse  images  in  (10). 

The  function  ri(y)  is  obtained  by  RofL(Y);  RofL(Y)  calls  two  different  routines,  EC(Ind,Y) 
and  CP(Y).  The  function  EC(Ind,Y)  will  compute  E[ X  |  Yi  =  y\  shown  in  (19).  The  second 
part  of  (9),  which  depends  on  the  {c*,},  is  computed  by  CP(Y).  To  compute  P (V*  €  Akj  |  Yi  = 
y)  defined  in  (18).  SumcP(K,  J,Y)  is  used.  Since  a  given  subpartition  Akj  could  be  a  union  of 
intervals,  SumcP  computes  the  conditional  probability  at  each  interval  using  Fun4(LL,UL,Y) 
and  adds  the  results  together. 

Before  explaining  how  the  set  defined  in  (10)  was  obtained  to  best  characterize  the 
partition  A,  we  discuss  two  routines  that  were  used  in  this  process,  Levels  and  Root.  The 
routine  Levels  computes  the  different  levels  needed  to  generate  the  partition  f°r  a 

given  sensor  k.  If  there  are  N  subpartitions  we  will  have  N  —  1  levels  all  stored  in  an  array 
Lev(I)  defined  as: 


Lev(I) 


C(sensorindex,  I) +  C(  sensor index,  I+l) 

2 


where  I  —  1, . . . ,  jV  —  1. 

Given  the  function  RofL(Y)  and  an  interval  [LL,UL]  the  routine  Root  will  find  a  simple 
root  in  this  interval.  If  more  than  one  root  exists  in  a  given  interval,  an  error  message  is 
issued.  To  compute  the  roots  the  NAG  routine  C05ADF  is  used. 

To  obtain  a  new  partition  A,  two  main  routines  are  used:  ListofZero  and  Nevpart. 
Listof  Zero  will  find  all  the  zeros  of  Rof  L(Y)  —  Lev(I).  To  do  this,  we  sample  r/(y)  for  Numpt 
values  of  y.  The  values  of  y  are  stored  in  the  array  Yp.  The  values  of  r/(y)  are  stored  in 
the  array  Rp.  The  routine  ListofZero  will  search  for  a  change  of  sign  of  Rp(J)-Lev(I). 
A  change  of  sign  is  assumed  to  indicate  that  at  most  one  root  exists  in  the  interval  Y(J) 
and  Y(J+1).  This  assumes  that  enough  samples  of  r/(y)  are  taken.  The  zeros  are  stored 
in  the  array  Zero(I,K),  where  the  index  I  indicates  to  which  level  the  zero  corresponds. 
In  addition,  the  nearest  value  to  the  left  and  to  the  right  of  r/(y)  are  stored  respectively 


c. 


in  Testpt(I,K,l)  and  Testpt(I ,K,2).  The  number  of  zeros  for  each  level  is  stored  in  the 
array  NumofZero. 

The  routine  Newpart  will  use  the  list  of  the  zeros  of  the  different  levels  and  the  test  points 
to  obtain  the  new  partition  A.  The  test  points  will  be  used  to  take  a  peak  at  the  left  or 
right  of  a  zero  and  will  enable  us  to  determine  if  the  range  of  r/(y)  in  an  interval  formed  by 
two  consecutive  zeros  lies  above  or  below  a  certain  level  or  between  two  levels.  To  determine 
the  intervals  of  y  for  which  r/(y )  is  between  two  levels,  we  will  first  merge  in  order  using 
MergeSort  the  zeros  and  their  corresponding  test  points  of  the  two  levels  in  the  arrays  Mzero 
and  Mtestpoint  respectively  and  then  find  the  desired  intervals. 

Finally  the  routine  UL(k,AR)  where  k  =  l,...,n  will  return  a  new  partition  of 
denoted  by  {A*,}^;  and  stored  in  the  array  AR  given  the  N  x  n  array  C  that  contains  the 
different  elements  c*.,.  This  routine  is  used  in  the  main  program.  It  plays  the  role  of  Ui(C ,  A) 
defined  in  Definition  2  and  used  in  Algorithm  1. 

To  study  the  performance  of  Algorithm  1  we  computed  the  mean  square  error  E[  \X  —  X|2  ] 
defined  in  (11)  each  time  the  partition  A  and  its  corresponding  coefficients  c*.  are  updated. 
The  latest  mean  square  error  is  compared  with  the  previous  value.  In  case  of  tiny  substantia] 
improvement  another  iteration  is  taken,  otherwise  we  stop  running  the  program. 


Appendix  A.  A  Listing  of  the  Code 


c 


define. f 


c  This  variable  stores  the  Mean  Square  of  the  distribution 
c  of  X  U~(-A,A). 

DOUBLE  PRECISION  VARX 
COMMON/BLOOC45/VARX 


c  NS  --->  nunber  of  sensors 
c  N  --->  ixnfcer  of  partitions 

c  HAX*2*SUBP,  where  SUBP  is  the  nunber  of  subpartitions 
c  The  data  of  the  variables  NS, N, SUBP  are  read 
c  respect i vel ly  from  the  file  init.dat. 

INTEGER  NS, N, MAX 
COMMON/BLOCK5/NS , N , MAX 


c  A1  =NS*N  and  B1=NS*NS*N*N 
INTEGER  A1.B1 
COMMON/BLOCK6/A1 ,B1 


c  COND  — >  indicates  if  the  conditional  probabilities  are 
c  conditioned  on  Sensorl  or  Sensor2  (Z1  or  22). 

c  SNUMB — >  indicates  the  SENSOR  NUMBER  we  are  dealing  with 
c  If  we  want  the  LLOYD-MAX  partition  for  SENSOR  2 

c  then  SNUMB  is  2. 

c  If  we  are  running  the  algorithm  then  SNUMB=NS. 

INTEGER  COND, SNUMB 
COMMON/BLOCK7/COND 
COMMON/BLOCK25 /SNUMB 


c  NNS  --->  MAXIMUM  nunber  of  sensors 
c  NN  --->  MAXIMUM  nunber  of  partitions 
c  MSUBP -  -  - >  MAXIMUM  nurtaer  of  siiopertitions 
INTEGER  NNS , NN , NSUBP , NMAX , AS , NA 

PARAMETER  ( NNS=2 , NN=64 , NSUBP=30,NMAX=2*NSUBP , AS=NNS*NN , NA=AS*AS ) 


c  The  variable  LEVEL  is  needed  to  find  all  the  values  of 
c  y  for  which  the  equation  LEVEL=Rof L(y)  holds. 

DOUBLE  PRECISION  LEVEL 
COMMON/BLOCK23/LEVEL 


c  LMAX  -vindicates  if  we  are  computing  the  LLOYD-MAX  partition 
c  or  not, LMAX  =.TRUE.  if  we  are  ,LMAX». false  if  we  are  not 

LOGICAL  LMAX 
C0MM0N/BL0CK24/LMAX 


c  Rp  is  an  array  that  contains  the  values  of  RofL  at  specific 
c  values  of  y.Yp  is  an  array  that  contains  the  corresponding 
c  values  of  y  for  Rp. 

c  The  MAXIMUM  nunber  of  points  at  which  RofL  can  be  evaluated 
c  is  MAXPT. 

INTEGER  MAXPT 
PARAMETER  <MAXPT*6000) 

DOUBLE  PRECISION  Rp(MAXPT),Yp(MAXPT) 
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COMMON/BLOCK1 OO/Rp, Yp 


c  NUHPT  indicat**  th*  nuaber  of  point*  (Yp,Rp)  is  taken 
INTEGER  NUMPT 
CQMHON/BLOCK1Q1/NUMPT 


c  Each  level  has  a  certain  nunber  of  ZEROS  stored  in  the 
c  array  ZERO.ZN  Mill  be  the  maximun  nuiber  of  levels 
c  allowed  and  ZMAX  is  the  maximun  nuiber  of  zeros  for 
c  a  given  level. 

INTEGER  ZN.ZMAX, INTER 
P ARAMETER ( ZN=NN , ZMAX *30 , I NTER=2 ) 


INTEGER  HMAX, WINTER 
PARAMETER (MMAX=2*ZMAX ,M1NTERZ2) 


integer  azero,  nWcol,  nYcol 
parameter  (azero=0,  nWcol=5,  nYcol=2) 

c 

integer  lower,  upper,  meanl,  ms,  height 
parameter  (lower=1,  upper=2,  mean1=3,  ms=4,  height=5) 

c 

double  precision  Uinfo(azero:NNS,nWcol ),  Yinfo(1 :NNS, nYcol ) 
double  precision  111,  ull,  112,  ul2,  yy 
comnon  /select/  Uinfo,  Yinfo,  111,  ull,  112,  ul2,  yy 
c 

integer  igtype,  si,  s2 
common  /iselct/  igtype,  si,  s2 
c 

integer  corr,  jpr,  unipr 
parameter  (corr*1,  jpr=2,  unipr=3J 
c 

integer  cndprb,  cndexp,  convol 
parameter  (cndprb*1,  cndexp=2,  convol =3) 

c 

double  precision  relerr 
common/block26/relerr 

c  firspt  is  the  first  point  at  which  Rofl  is  evaluated 
c  lastpt  is  the  last  point  at  which  Rofl  is  evaluated 
double  precision  firstpt, lastpt 
common/block27/f irstpt, lastpt 

c  constant  MC 

double  precision  MC.varx.varz 
common/block28/NC,varz 
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c 


MAINU.F 


c  . MAIN . MAIN . MAIN . MAIN . MAIN 

c  This  is  the  Min  program. First  the  program 
c  calls  the  routine  INITDI,  this  routine  provides  the 
c  information  of  the  different  distributions  that  are 
c  needed  and  stored  in  arrays  available  to  other  routines, 
c  1)If  0  is  entered  the  algorithm  is  run  with  a  nutter  of 
c  NIT  iterations, if  we  want  to  continue  where  we  left  off 
c  (after  NIT  iterations), the  last  partition  generated  by 
c  the  program  is  made  available  in  file  "startl". 
c  2) If  1  or  2  is  entered  the  program  will  find  the  LLoyd-Max 
c  partition  for  sensor  1  or  2  those  partitions  are  available 
c  in  "startl"  or  "start2“. 
c  This  is  the  following  steps  taken  by  MAIN  : 
c  1)  The  initial  partition  is  stored  in  the  array  AR 

c  using  ReadAR(AR). 

c  2)  The  routine  SN(AR,C,MSE)  is  used  to  compute  the 

c  optimal  coefficients  for  a  given  partition 

c  and  are  stored  in  the  array  C.This  routine 

c  also  computes  the  mean  square  error  stored  in  MSE. 

c  3)  To  find  a  better  partition  we  iterate  a 

c  a  fixed  nutter  of  times  (NIT). After  each 

c  iteration  an  updated  partition  along  with  its 

c  optimal  coefficients  is  computed  using  UL(COND,AR) 

c  and  SN(AR,C,MSE)  respectively, 

c 

c  The  program  will  used  two  input  files. 

c  1)  distr.dat  — >  Contains  the  distribution  informations 

c  2)  init.dat  --->  Indicates  NS, N, SUSP, REIERR, NIT, TEST  and 

c  DELTA  respectively. 

PROGRAM  main 

INCLUDE  'define. f' 

INCLUDE  ' Camay. f' 

INCLUDE  'ARarray.f ' 

INCLUDE  'MEANarray.f ' 

INCLUDE  'NUarray.f ' 

INCLUDE  'TEMP2array.f ' 

INTEGER  COUNT, TEST, I, J,K, IA, NIT, NUMBPT,NPT,KL 
DOUBLE  PRECISION  Yx,DUM,X1 ,X2,Y1 ,Y2, SUSP, MSE 
DOUBLE  PRECISION  MSI .DELTA, RofL.PT 
EXTERNAL  ClearAR.ReadAR.Ul , MSI .LOAD, SN, INITDI 
EXTERNAL  CreateRofL.RofL 


OPEN  (UNIT=81, FILE® 'init.dat '.STATUS® 'old') 

READ(81 ,*)  NS, N, SUBP, REIERR 
REA0(81 ,*)  NIT, TEST, DELTA 
READ(81 ,*)  MC 
CLOSE  (UNIT=81 ) 

WRITE(*,*)  'Please  Enter  the  Nutter  of  ITERATIONS 
WRITE(*,*)  'NIT®',  NIT 
c  READ(*,*)  NIT 
URITE(*,») 

MAX®2*SUBP 

URITE(*,*)  'The  Maximua  nutter  of  Sensors  is  ',NNS 
WRITE(*,*)  'The  Maximua  nuiter  of  Partitions  is  ' ,NN 
URITE(*,*)  'The  Maxim*  nuiter  of  Subpartitions  is  '.NSUBP 
WRITE(*,*)  'The  Nuiter  of  Sensors  is  ',NS 
WRITE(*,*)  'The  Nuiter  of  Partitions  is  ',N 
WRITE(*,*) 
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CALL  INITOI 


WRITE!*,*) 

WRITE!*, *)'ENT£R  0  to  RUN  the  ALGORITHM' 

WRITE!*, *)'ENTER  1  to  obtain  the  LLOYD-MAX  for  SENSOR  V 
WRITE!*,*) 'ENTER  2  to  obtain  the  LLOYD-MAX  for  SENSOR  2' 
WRITE!*,*)  'TEST-',  TEST 
c  READ!*,*)  TEST 
X1=Winfo(ZERO,MS) 

X2=Winfo(TEST,MS) 

Y1=Winfo!ZERO,MEAN1) 

Y2=Winfo!TEST,MEAN1 ) 

IF  (TEST.EO.O)  THEN 
LMAX=. FALSE. 

VARX-X1 

ELSE 

LMAX=.TRUE. 

COND=TEST 

NS=1 

VARX=X1  ♦  X2  ♦  2*Y1*Y2 

END  IF 

A1=NS*N 

B1=A1«A1 

CALL  ClearAR(AR) 

CALL  ReadAR  (AR) 


WRITE!*,*)  '  The  RELERR  is: ' ,RELERR 
CALL  SN(AR,C,MSE) 

WRITE!*,*)  'The  initial  MEAN  SQUARE  ERROR  is:  ',MSE 
WRITE!*,*) 

DO  11  11=1, NS 

If  ( LMAX )  THEN 
XL=T£ST 
ELSE 
KL=K 
END  IF 

X 1 =Y inf of KL, LOWER) 

X2=Y inf of XL, UPPER) 

WRITE!*, 14)  XL 
WRITE!*, 12)  XL,X1 ,X2 

11  CONTINUE 
WRITE!*, 13)  DELTA 

14  FORMAT! 'For  Sensor' , 12, ':' ) 

12  FORMAT ! 3X , 1 Rof ' , 1 1 , ' ! Y )  sampled  for:  ',F7.4,'  <  Y  «',F7.4) 

13  FORMAT! 'The  INCREMENT  DELTA  used  for  both  cases  is:',F8.5,/) 
COUNT  =0 

DO  55  Ia=1 ,NIT 
COUNT =COUNT+ 1 

WRITE!*,*)  '********* '.count,'  ITERATION***************' 

DO  15  X=1 ,NS 

IF  fLMAX)  THEN 
CONO-TEST 
SNUMB-1 
ELSE 

COND-k 

SNUMB-k 

END  IF 
J=1 

X 1 =Y inf of CONO, LOWER) 

X2=Yinfo(COND, UPPER) 

PT=!X2-X1)/P£LTA 
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NUMBPT*OINT(PT) 

NPT=NUMBPT  +1 

OPEN ( UN  I T=2 1 , F I LE* ' out 1 ' , STATUS* ' UNKNOWN ' ) 
c  WR1TE(*,*)  'Oat*  for  Rofl  is  being  computed..' 

00  10  1*1, NPT 

Yx=X1  ♦  DELTA* ( 1 -1 ) 

IF  ((Yx.GT.X1) .ANO. (Yx.LT ,X2) )  THEN 
Rp( J)*RofL(Yx) 

Yp( J)*YX 

IF  (J.LE.MAXPT)  THEN 
WRITE(21 ,32)  Yx,Rp(J), J, I 

32  F0RNAT(F13.6,2x,F13.6,2x, 14, 2x, 14) 

ELSE 

WRITE(*,*)  'ERROR  IN  DIMENSION  FOR  Rp  ARRAY' 
ENO  IF 

j=j+1 

ENO  IF 

10  CONTINUE 

c  wri te(*,*)  'ENTER  ANY  NUMBER  TO  /NTINUE . ' 
c  read(*,*)  OUM 

CLOSE  (UNIT=21 ) 

NUMPT*J- 1 

CALL  Ul( SNUMB , AR ) 

15  CONTINUE 


CALL  SN(AR,C,MSE) 

URITE(*,*)  'The  MEAN  SQUARE  ERROR  is:  ',MSE 
WRITE(*,») 

55  CONTINUE 

STOP 

END 
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GEN.  F 


c  _ read  data  of  partition  into  a  3-D  arrray  AR.. 

c  This  routine  read  fro*  a  file  "start*1  the  data  of  a  given 
c  partition  and  store  into  a  3-D  array  AR.The  index  I  of  the 
c  matrix  AR  represent  the  sensor  nuifcer.the  index  J  represent  the 
c  the  subpartition  nunber  for  a  given  sensor  1  finally  the  indices 
c  K  and  K+1  indicate  the  lower  and  upper  values  respectively  of  an 
c  interval  of  the  siipartition  J  of  a  sensor  l.A  subpartition  J  of 
c  a  given  sensor  I  can  be  a  union  of  intervals. When  reading  from 
c  "start"  each  subpartition  J  is  separated  by  the  two  consecutive 
c  nunbers  of  - 1 . 

SUBROUTINE  ReadAR  (AR  ) 

INCLUDE  'define. f' 

DOUBLE  PRECISION  AR(NNS,NN,NHAX),LL,UL 
INTEGER  I,J,K 

OPEN  (UNI  T  =  10, FILE='start ' , STATUS='old' ) 

DO  10  1=1, NS 

DO  20  J=1,N 
K  =  1 

30  CONTINUE 

READOO,*)  LL ,UL 
AR( I , J , K)=LL 
AR( I , J ,K*1 )=UL 
K=K+2 

IF  ( (UL .EO. - 1 . 000) .AND . (LL .EO. - 1 .000) )  THEN 
AR(I,J,K-2)=0.D0 
AR( I , J ,K- 1 )=0 .DO 
GO  TO  20 
END  IF 
GO  TO  30 

20  CONTINUE 

10  CONTINUE 

RETURN 
END 


c  . INITIALIZE  AR  TO  ZERO 

SUBROUTINE  ClearAR  (AR) 

INCLUDE  'def ine.f ' 

INTEGER  !,J,K 

DOUBLE  PRECISION  AR ( NNS , NN , NMAX ) 

DO  100  1=1, NS 

DO  200  J=1 , N 

DO  300  K=1 ,MAX 
AR( I ,J,K) =0.000 
300  CONTINUE 

200  CONTINUE 

100  CONTINUE 
RETURN 
END 


C  This  routine  returns  the  matrix  NU(I,J).For  example  if 

c  N(I,J)=10  it  means  that  a  subpartition  J  of  a  given  sensor 
c  I  is  formed  from  the  union  of  five  different  intervals. 


c  This  matrix  wilt  be  used 
SUBROUTINE  COUNTP  (NU) 

INCLUDE  'define.#' 
INCLUDE  'ARarray.f ' 

INTEGER  NU(NHS,NH) 

INTEGER  L.J.KB 


DO  20  L=1 , NS 
DO  30  J=1,N 
NUCL, J)=0 
KB=1 

10  CONTINUE 

IF  (AR(L,J,KB).NE.AR(L,J,KB+1))  THEN 
IF  (KB.LT.MAX)  THEN 
NU(L,J)=NU(L,J)  +2 
KB=K8+2 
GO  TO  10 

END  IF 
END  IF 

30  CONTINUE 

20  CONTINUE 
RETURN 
END 


c  Given  two  subpartitions  Aki  and  A l j  Eli  (Yk  )I  <Y l >]  is  computed 
c  as  follows: 


c 

c 

c 

c 

c 

c 

c 

c 

c 


1-  If  k=l  and  i=j  then 

PlYk  Aki)  is  computed. Si  nee  Aki  the  subpartition 
I  of  a  sensor  K  can  be  a  union  of  intervals  the 
probability  of  Yk  along  each  interval  is  computed 
and  then  sunned  together  to  perform  this  correctly 
the  nunber  of  intervals  was  needed  (NU). 

P(Yk  Aki)  was  computed  using  PG. 

2-  If  k=l  and  i  NE  j  then 

the  routine  is  assigned  the  value  ZERO 


c  3-  If  k  HE  l  then 

c  P(Yk  Aki.Yl  Alj)  is  computed, 

c  P(Yk  Aki.Yl  Alj)  was  computed  using  JPG. 

c  E tl  (Yk  )1  (Yl)]  is  computed  and  stored  in  the  matrix 

c  SUMP(L , J,K, I ) ,  this  matrix  is  needed  to  compute  the  coefficients 
c  contained  in  the  matrix  C. 


c  COMPUTE  SUM  OF  PROBABILITIES  IF  WE  HAVE  UNION  OF  INTERVALS... 

OOU8LE  PRECISION  FUNCTION  SUMP(L , J ,K, I ) 

INCLUDE  'def ine.f ' 

INCLUDE  'ARarray.f' 

INCLUDE  'NUarray.f ' 

DOUBLE  PRECISION  RES.ANS 
INTEGER  NU 1 , NU2 , NX , MX , I , J , K , L 
DOUBLE  PRECISION  JPG.PG 
EXTERNAL  JPG.PG 

NU1=NU(L, J) 

NU2=NU(K, I ) 

RES=0.000 
IF  (L.EQ.K)  THEN 
IF  (I.EO.J)  THEN 

DO  15  NX*1,NU1,2 

ANS=PG(AR(L, J,NX),AR(L,  J.NX+1 ), L) 
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RES=RES+ANS 
IS  CONTINUE 

SUMP=RES 

ELSE 

SUNP=0.00 
END  IF 
ELSE 

DO  30  MXS1 ,NU1 ,2 

oo  40  nx»i,nu2,2 

ANS=JPG(AR(L, J,MX),AR(L, J.HX+I ) ,AR(X, I ,NX) , 
*  AR(K, 1 ,NX+1 ),L,K) 

R£S=RES+ANS 
40  CONTINUE 

30  CONTINUE 

SUMP=RES 

END  IF 
RETURN 
END 


c  This  routine  will  compute  E[XI  (Yk)]  the  integral  of 
c  fyk(y)*E lX/Yk=Y]  over  the  subpartition  Aki.lt  will  use 
c  MEANX  a  routine  that  computes  the  integral  of  fyk(y)*E CX/YkI 
c  over  a  given  interval.  Note  that  Aki  is  formed  of  a  union  of 
c  intervals  .knowing  those  intervals  (AR).their  number  (NU)  and 
c  using  MEANX  E[XI  (Yk)]  is  computed  and  stored  in  the  matrix 
c  SUMM(I.J)  index  I  indicating  the  sensor  and  J  the  subpartition 
c  SUMM  is  needed  to  compute  the  coefficients  contained  in  the 
c  matrix  C. 
c 

c  . sum  of  means . 

DOUBLE  PRECISION  FUNCTION  SUMM(I.J) 

INCLUOE  'def ine.f ' 

INCLUDE  'ARarray. f ' 

INCLUOE  'NUarray. f ' 

INTEGER  NU1 , NX, I , J 
DOUBLE  PRECISION  RES, ANS 
DOUBLE  PRECISION  MEANX 
EXTERNAL  MEANX 

NU1=NU(I ,J) 

RES=0.000 

DO  15  NX=1 , NU1 , 2 

ANS=MEANX(AR( I , J ,NX) , ARC  I , J ,NX+1 ), I ) 

RES=RES+ANS 
15  CONTINUE 
SUMM=RES 
RETURN 
END 


c  . generate  the  R  matrix . 

c  This  routine  will  generate  the  R  matrix  (a  symmetric  matrix) 
c  needed  to  solve  the  optimal  coefficients, 
c  The  equation  that  is  needed  to  be  solved  is: 
c  Rc=mean 

c  where  c  is  the  matrix  that  contains  the  optimal  coefficients. 
SUBROUTINE  GENR  (R) 

INCLUOE  'define. f' 

INCLUDE  'ARarray. f' 

INCLUOE  'NUarray. f' 

INCLUOE  'TEMP2array.f ' 
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INTEGER  II.JI.I.J.IC.L 
DOUBLE  PRECISION  R(AS,AS) 
OOUBLE  PRECISION  SUMP 
EXTERNAL  SUMP 


11=0 

DO  10  L=1,NS 

DO  20  J=1,N 
11=11+1 
J1=0 

DO  30  K=1 ,NS 
DO  40  1=1 ,N 

J1=J1+1 

IF  (J1.GE.I1)  THEN 

R ( 1 1 , J 1 ) =SUMP ( L , J , K, I ) 

IF  (R( 1 1 , J1 ) .NE .O.OdO)  THEN 
c  wri te(*, 100)  R(11,Jl),I1,Jl 

c  100  FORMAT  (F16.12,  '  R( ' , 13, ' , ' , 13, ' ) ' 5 

END  IF 

ELSE 

R(  II , J1 )=R( J1 , 1 1 ) 

END  IF 

40  CONTINUE 

30  CONTINUE 

20  CONTINUE 

10  CONTINUE 
RETURN 
END 


c  . generate  mean  matrix . 

c  This  routine  will  generate  the  other  matrix  needed 
c  to  solve  for  the  optimal  coefficients  the  mear  natrix. 
c  The  equation  that  is  needed  to  be  solved  is: 
c  Rc=mean 

c  where  c  is  the  matrix  that  contains  the  optimal 
c  coefficients. 

SUBROUTINE  GENM  (MEAN) 

INCLUDE  'define. f' 

INCLUDE  'ARarray.f ' 

INCLUDE  'NUarray.f ' 

INTEGER  IND,J,L 

DOUBLE  PRECISION  MEAN(AS) 

DOUBLE  PRECISION  SUMM 
EXTERNAL  SUMM 

I ND = 1 

DO  10  L=1,NS 

DO  20  J=1,N 
MEAN(IN0)=SUMH  (L,J) 

c  wri te(*, 100)  MEAN(IND), INO, l, j 

c  100  format(F16.6, '  MEAN( ' , 13, ' ) ' ,213) 

IND=IND+1 
20  CONTINUE 

10  CONTINUE 
URITE(*,+) 

RETURN 

END 


c  . SOLVE  LINEAR  EOT . 

c  This  routine  will  solve  the  equation:  Rc=Mean 
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c  c  Mill  be  the  matrix  that  contains  the  optimal 
c  coefficient. To  obtain  those  coefficients  the 
c  SVD  approach  was  used, 
c  The  NAG  routine  F02WEF  used  performs  on  R 
c  the  matrix  factorization  known  as  the  singular  value 
c  decomposit ion, SOLVE  also  computes  the  minimized  mean 
c  s«Mare  error  stored  in  MSE. 


SUBROUTINE  SOLVE  ( R, MEAN, C, MSE) 

INCLUDE  'define. f' 

INTEGER  LOR , LDM, LDPT ( NCOLM 

PARAMETER  ( LDR=AS , LDM=AS , LDPT=AS , NCOLM* 1 ) 

INTEGER  LWORK 

PARAMETER(LW0RK=AS**2  +  5*(AS- 1 ) ) 

INTEGER  I.IFAIL, J.IND 
LOGICAL  UANTP.WANTQ 


DOUBLE  PRECISION  R(LDR,AS),MEAN(AS) ,C(NNS,NN) ,C1 (AS) , DUMMY (1 ) 
DOUBLE  PRECISION  UORK( LWORK ) , SV( AS ) , PT ( LDPT , AS ) , MEANT ( LDM) 
DOUBLE  PRECISION  RES, MSE, TAU(S, F06EAF 
EXTERNAL  F06EAF , F02WEF 

DO  30  1*1, AT 

MEANT ( I )=MEAN( I ) 

30  CONTINUE 
UNIQUE*. TRUE. 

WANTQ=.TRUE. 

WANTP=.TRUE. 

I  FA  I L  =0 

CALL  F02WEF  (A1 ,A1 ,R, LDR, NCOLM, MEANT ,L0M,WANTQ,0UMMY , 1 ,SV, 

*  UANTP.PT, LDPT, WORK, I  FAIL) 

IF  ( I FAIL.NE.O)  THEN 

wri te(*,*)  '  ERROR  IN  SVD  DECOMPOSITION' 

END  IF 

OPEN  (UN I T= 1 , F I LE= ' svd . dat ' , STATUS* ' unknown ' ) 

WRITE(1 , 1 1 ) 

WR I TE( 1 , 101 )  (SV( I ), 1=1 ,A1 ) 

URITE( 1 , 12) 

DO  10  1=1, A1 

c  WRITE(1 , 101 )  (PT( J, I ) , J=1 , A1 ) 

10  CONTINUE 
WR  I TE ( 1 , 13) 

DO  20  1=1 ,A1 

c  URITE( 1,101)  (R(I,J),J*1,A1) 

20  CONTINUE 

URITE(1 , 14) 

WRITE(1 , 101 )  (MEANT ( I ) , 1=1 ,A1 ) 

CLOSE  (UNIT=1) 

c  TAU  is  the  absolute  error  tolerance 
TAU=RELERR*SV(A1 ) 

DO  65  1=1, A1 

IF  (SV(l).LT.TAU)  GOTO  75 
NSV=I 

65  CONTINUE 

75  CONTINUE 

write(*,*)  NSV,'  singular  values  are  >=',TAU 
Hri te(*,*)  'SVMIN*  ',SV(1) 
write(*,*)  'SVLAST*  ',SV(NSV) 
write(*,*)  'SVMAX*  <,SV(A1) 

DO  50  1=1 ,A1 

C1(I)=0.000 
DO  40  J=1 ,NSV 


3S 


S=M£ANT (J)/SV(J) 

Cl C I )=C1< I )  *PT(J,I)*S 
40  CONTINUE 

SO  CONTINUE 

c  WRITE!*,*)  'The  solution  to  the  system  is:' 
c  WRITE!*, 100)  (Cl ! I ) , 1=1 ,A1 ) 


c  Dot  Product  of  the  two  vectors  Cl  and  MEAN 
RES=F06£AF(AS,C1 , 1 .MEAN , 1 ) 

MS£=VARX-RES 


IN0=1 

00  15  1=1, NS 

DO  25  J=1,N 
C(  I , J  )=C1 ( INO) 

IND=IND+1 
25  CONTINUE 

15  CONTINUE 

100  FORMAT(4(1X,D12.4)) 

11  FORMAT!/'  SINGULAR  VALUE') 

12  FORMAT!/'  RIGHT-HANO  SINGULAR  VECTORS  BY  COLUMN') 

13  FORMAT!/'  LEFT-HANO  SINGULAR  VECTORS  BY  COLUMN') 

14  FORMAT!/'  VECTOR  0"*MEAN  ') 

101  FORMAT(5( 1X.D12.4) ) 

RETURN 

END 


c  This  routine  generates  the  matrix  c  that  contains 

c  the  optimal  coefficients  given  a  partition  A. 
c  This  routine  is  used  in  the  MAIN  program. 

SUBROUTINE  SN(AR,C,MSE) 


INCLUDE  'define. f' 

INCLUDE  'MEANarray.f ' 

INCLUDE  'NUarray.f ' 

DOUBLE  PRECISION  R(AS, AS) , C(NNS,NN) , AR(NNS,NN,NMAX) ,MSE 
EXTERNAL  COUNTP , GENM ,  GENR ,  SOLVE 

CALL  COUNTP(NU) 

CALL  GENM(MEAN) 

CALL  GENR(R) 

CALL  SOL VE (R, MEAN, C.MSE) 

RETURN 

END 


c  . sum  of  conditional  prob . 

c  This  routine  is  needed  to  compute  Rl(Y). 
c  It  computes  the  conditional  probabilities  for 
c  a  given  subpartition. 

DOUBLE  PRECISION  FUNCTION  SUMCP(I,J,Y) 

INCLUDE  'define. f' 

INCLUDE  'ARarray.f' 

INCLUDE  'NUarray.f' 

INTEGER  I , J, IND 
DOUBLE  PRECISION  RES.ANS, Y 
INTEGER  NU1.NX 
DOUBLE  PRECISION  FUN4 


EXTERNAL  FUN4 


NU1=NU( I , J) 

RESZ0.000 

INDzCOND 

DO  15  NX>1 ,NU1 ,2 

ANS=FUN4(AR( I , J,NX) ,AR( I , J ,NX+1 ) , Y,  I , IND) 
c  write<*,*>  'arf'.i.'.'.j.'.'.nx,')',  ar(i,j,nx> 

c  write(*,*)  'ar<',i,V.i.V,nx+V>'.  arc  t .  j  ,rvx+1  > 

RES*RES+ANS 
15  CONTINUE 
SUMCP ZRES 

c  Write(*,*)  res,'  res  ' 

RETURN 

ENO 


c  This  routine  is  needed  to  compute  Rl(Y). 

DOUBLE  PRECISION  FUNCTION  CP  (Y) 

INCLUDE  'define. f' 

INCLUDE  'Camay. f' 

INCLUDE  'ARarray.f ' 

INCLUDE  'NUarray.f ' 


INTEGER  I ND , J ,  K 
DOUBLE  PRECISION  Y 
DOUBLE  PRECISION  RES.ANS 
DOUBLE  PRECISION  SUMCP 
EXTERNAL  SUMCP 

IND-SNUMB 
RES=0.0D0 
DO  10  K=1,NS 

IF  (K.NE.IND)  THEN 

DO  20  J=1,N 

ANS=SUMCP  (K, J,Y)*C(K, J) 

RES=ANS+RES 

c  write<*,*)  c(k, j ),  sumcp<k, j,y),'  c<k,j>' 

20  CONTINUE 

END  IF 

10  CONTINUE 
CP=RES 

c  write(*,*)  res,'  fres' 

RETURN 

END 


c  .  Rl(Y) . . . 

c  This  routine  computes  Rl(y). 

DOUBLE  PRECISION  FUNCTION  RofL(Y) 

INCLUDE  'define. f' 

INCLUOE  'ARarray.f' 

INCLUDE  ' Camay. f' 

INCLUDE  'NUarray.f' 


INTEGER  I  NO 
DOUBLE  PRECISION  Y 
DOUBLE  PRECISION  CP, EC 
EXTERNAL  CP, EC 


I ND=COND 

c  IF  (LMAX)  THEM 
Rof L=Y 

c  ELSE 

RofL=EC( IND, Y)  -CP(Y) 
c  ENO  IF 
RETURN 
ENO 


c . RofLlev . 

DOUBLE  PRECISION  FUNCTION  RofLlev(Y) 

INCLUDE  'define. V 

DOUBLE  PRECISION  Y.RofL 
EXTERNAL  Rof L 

RofL lev=RofL(Y)  -  LEVEL 

RETURN 

END 


c 


ROOTFB.F 


c  This  function  finds  one  of  the  sinple  roots  of  R(Y)  given 
c  the  interval  [ll,ul].  If  more  than  one  root  exists  in  the 
c  given  interval  an  error  message  is  issued, 

c  The  subroutine  C05ADF  from  the  HAG  library  is  used. 

DOUBLE  PRECISION  FUNCTION  ROOT(LL,UL,RofL) 

DOUBLE  PRECISION  LL,UL, RofL.EPS, ETA, X 
EXTERNAL  RofL,C05ADF 
INTEGER  [FAIL 

EPS=1 .00 -8 

ETA=0.000 

IFA1L=1 

CALL  C05ADF(LL,UL,EPS,ETA,Rof L,X, I  FAIL) 

IF  (IFAIL.GT.0)  THEN 

WRITE<*,*)  I  FA I L  , '  ERROR  CHECK  C05DAF  in  ROOT' 

END  IF 
ROOT=X 
RETURN 
END 


c  This  routine  computes  the  different  levels. 

c  The  assunption  in  this  algorithm  is  that  C(l,1)<. _ <CC l , N > 

c  this  is  why  the  different  elements  of  the  n  X  N  matrix  C 
c  are  sorted  in  a  increasing  order  before  the  computation  of 
c  the  corresponding  levels. 

SUBROUTINE  LEVELS  (LEV) 

INCLUDE  'define. f' 

INCLUDE  'Carray.f' 

INTEGER  DUM,K,J,IND,I 

DOUBLE  PRECISION  LEV(NN),C(NNS,NN),TEHP 


0UM=N- 1 
IND=SNUMB 
DO  15  K=1,N-1 

DO  25  J=N,K+1 , -1 

IF  (C(IND,J).LT.C(IND,J-1))  THEN 
TEHP=C( IND , J  - 1 ) 
C(IND,J-1)=C(IND,J) 
C(IND,J)=TEMP 

END  IF 

25  CONTINUE 

15  CONTINUE 

DO  10  1=1, DUM 

LEV( I)=(C(IND,l*1 )+C( IND, I )}/2.0DO 
10  CONTINUE 
c  write(*,*) 

write<*,*)  'The  levels  used  are:  ' 
urite(*,100)  (lev(i),i=1,dua) 

100  format(4F16.6) 

uritet*,*) 

RETURN 

END 


c  This  routine  coa^utes  all  the  zeros  of  Rp(Y)  at  a  given 
c  level  lev(!).The  zeros  are  stored  in  the  matrix  ZERO  also  the 

c  nearest  value  of  Rp(Y)  to  the  left  and  to  the  right  of  t 
c  the  zero  is  stored  in  the  matrix  TESTPT, finally  the 


c  number  of  zeros  are  stored  in  NUMOFZERO. 
c  To  find  the  zeros  the  routine  scans  a  list  of  data  that 
c  contains  values  of  Rp(Y)-levO)  in  search  of  a  change  of 
c  sign  or  a  zero  value. 

SUBROUTINE  USTOFZERO  ( ZERO, NUHOFZERO, TESTPT, LEV) 

INCLUDE  'define. f* 

INTEGER  I , J ,K,HAXPT 
INTEGER  NUMOFZERO(ZN) 

DOUBLE  PRECISION  TESTPT(ZN,ZMAX, INTER), LEV(NN), ZERO  (ZN.ZKAX) 

DOUBLE  PRECISION  ROOT, CHECK 

DOUBLE  PRECISION  RofLlev 

EXTERNAL  LEVELS, ROOT 

EXTERNAL  RofLlev 


DO  10  1=1, N-1 
K=1 

LEVEL=LEV( 1 ) 

DO  20  J=1 , NUMPT- 1 

CHECK=(Rp( J )-LEV( I ) )*(Rp( J+1 ) -LEV( I ) ) 

IF  (CHECK. LT.O.ODO)  THEN 

ZERO( I, K)=ROOT(Yp( J),Yp< J+1), RofLlev) 
TESTPT< I ,K, 1 )=Rp( J) 

TESTPT (1 ,K,2)=Rp(J+1) 

K=K+1 

ELSE  IF  (Rp(J).EO.LEV(I))  THEN 
IF  (J.NE.1)  THEN 
GO  TO  20 

ELSE 

ZERO( 1 ,K)=Yp(  1 ) 

TESTPT ( I ,K, 1 )=LEV( I )-Rp(  J+1 ) 
TESTPT ( I ,K,2)=Rp( J+1 ) 

K=K+1 
END  IF 

ELSE  IF  (Rp( J+1 ) .EO.LEVC I ) )  THEN 
ZERO( 1 ,K)=Yp( J+1 ) 

TESTPT ( I ,K, 1 )=Rp( J ) 

TESTPT ( I ,K,2)=Rp( J+2) 

K=K+1 

END  IF 

20  CONTINUE 

NUMOFZEROC I )=K* 1 
10  CONTINUE 
RETURN 
END 


c  This  routine  will  merge  in  order  two  lists  of  zeros  and 
c  testpoints.The  list  of  zeros  obtained  at  the  level  LEVA  is 
c  merged  with  the  list  of  zeros  obtained  in  LEVB  in  an  increasing 
c  order, the  new  list  of  zeros  is  stored  in  the  array  NZERO. 
c  Similarly  the  testpoints  obtained  at  the  two  different  levels 
c  are  stored  in  increasing  order  in  MTESTPT. 

SUBROUTINE  HERGEANDSORT  (LEVA, LEVB, ZERO, TESTPT, HZERO.NTESTPT, 

*  NUHOFZERO) 

INCLXE  'define. f' 

DOUBLE  PRECISION  ZERO(2N,2HAX) .TESTPT (ZN.ZKAX, INTER) 

DOUBLE  PRECISION  NZERO(NHAX) .NTESTPT (UMAX, WINTER) 

INTEGER  NUMOFZERO(ZN) 


INTEGER  LEVA, LEVS, ILA, ILB, IN, T 


ILA»1 

ILB=1 

IM*1 

10  CONTINUE 

IF  ( I LA . LE . NUNOFZEROt  LEVA ) )  THEN 
IF  (lLB.LE.NUMOFZERO(LEVB))  THEN 

IF  (ZEROCLEVA, ILA).LE.ZERO(LEV8, 1  LB ) )  THEN 

MZERO( IM)=ZERO(LEVA, ILA) 

MTESTPT (IN, 1 )=TESTPT(LEVA, I  LA, 1 ) 
MTESTPT  < IH,2)=TESTPT(LEVA, I  LA, 2) 
ILA=ILA+1 


ELSE 


MZERO( IM)=ZERO(LEVB, ILB) 

MTESTPT ( IN, 1 )=TESTPT(LEVB, ILB, 1 ) 
MTESTPT (IM, 2 )=TESTPT(LEVB, ILB, 2) 
ILB=ILB+1 


END  IF 
I M= I M+ 1 


GO  TO  10 
ENO  IF 
END  IF 

IF  ( ILA.GT.NUMOFZERO(LEVA))  THEN 
DO  20  T=ILB,NUMOFZERO(LEVB) 

MZEROC IM)=ZERO(LEVB,T) 

MTESTPT ( IM, 1 )=TESTPT(L£VB,T, 1 ) 
MTESTPT ( IM,2)=TESTPT(LEVB, T ,2) 
IM=IM+1 

20  CONTINUE 

ELSE 

DO  30  T=ILA,NUMOFZERO(LEVA) 

MZEROC IM)=ZERO(LEVA,T) 
MTESTPT(IM,1)=TESTPT(LEVA,T,1) 
MTESTPT ( IM, 2)=TESTPT(LEVA, T, 2) 
IM=IM+1 

30  CONTINUE 

END  IF 

RETURN 


ENO 


c  This  routine  given  a  list  of  the  zeros  and  its  testpoints 
c  obtained  at  a  given  level  gives  the  new  partition  stored  in 
c  the  matrix  ARS.The  testpoints  will  be  used  to  take  a  "peak" 
c  at  the  left  or  right  of  a  zero  which  will  enable  us  to 
c  determine  if  the  range  of  r(y)  in  an  interval  formed  by  two 
c  consecutive  zeros  is  above  or  below  a  certain  level 
c  or  between  two  levels. 

SUBROUTINE  NEUPART  (ZERO, TESTPT, INOX, NUNOFZERO, ARS, LEV) 

INCLUDE  'define. V 

INTEGER  !NOX,NUK,L,LEVA,LEVB,l,J,K 

DOUBLE  PRECISION  ZERO< ZN.ZMAX), TESTPT (ZN.ZMAX, INTER) 

INTEGER  NUHOFZERO(ZN) 

DOUBLE  PRECISION  ARS(NNS,NN,NNAX),LEV(NN) 

DOUBLE  PRECISION  HZERO(MMAX ) .MTESTPT (MMAX,MI NTER ) 

EXTERNAL  LEVELS, MERGEANOSORT 


c 

c 


L=INDX 
00  10  1=1, N 
K=1 

We  are  searching  for  all  the  values  of  Y  that  belongs  to 
[-B-A,B*A]  such  that  Rl<Y)  is  less  than  IEV(1). 

IF  CI.E0.1)  THEN 

00  20  J=1,NUM0FZER0(1) 

IF  (TESTPT ( I , J, 1 ) .LT.IEV(1 ))  THEM 
IF  CJ.E0.1)  THEN 

ARSCL, I ,IC)=  ZER0< I ,  J)+2 

ARSCL, I  ,IC+1 )=ZER0( I , J) 

writeC*,100)  ars(l, i ,k),arsCl, i , k+1 > , k 

K=K*2 

ELSE 

ARSCL, I ,K)=ZER0C I , J-1 ) 

ARSCL, I ,K+1 )=ZER0< I , J) 
wri teC*,100)  ars( l , i ,k),ars( l , i ,k+1),k 
K=K+2 
END  IF 
END  IF 

20  CONTINUE 

IF  C TESTPT C I , NUMOFZEROC 1 ),2).LT.LEV(1 ) )  THEN 
ARSCL, I ,K)=ZER0( I , J) 

ARSCL, l,K+1)=2EROCl,J)-3 

wri teC*, 100)  arsC l , i ,k),arsC l, i ,k+1),k 

K=K+2 

END  IF 

c  We  are  searching  for  all  the  values  of  Y  that  belongs  to 
c  C-S-A,B-A1  such  that  RICY)  is  greater  than  LEVCN-1). 

ELSE  IF  CI.EQ.N)  THEN 

00  30  J=1 , NUMOFZEROCN- 1 ) 

IF  CTESTPTCI*1,J,1).GT. LEVCN-1))  THEN 
IF  CJ.EQ.1)  THEN 

ARSCL,  l,K)=ZEROU-1,J)+2 
ARSCL, I ,K+1 )=  ZEROC I -1 , J) 
writeC*,100)  arsCl,i,k),arsCl,i,k+1),k 
IC=IC+2 

ELSE 

ARSCL, I ,K)=ZER0C I -1 , J- 1 ) 

ARSCL,  I  ,IC*1  )=ZEROC  I  -1 ,  J) 
writeC*,100)  arsCl,i,k),arsCl,i,k+1),k 
IC=IC+2 
END  IF 
END  IF 

30  CONTINUE 

IF  CTESTPTCl-1,NUMOFZERO(N-1), 2). GT. LEVCN-1))  THEN 
ARSC L , I ,K)=ZEROC I  - 1 , NUMOFZEROCN- 1 ) ) 

ARSCL, I ,K*1 )*ZEROC I  - 1 .NUMOFZEROCN- 1 ) )-3 
writeC*,100)  arsCl,i,k),arsCl,i,k+1),k 
Kaic+2 

ENO  IF 

c  We  are  searching  for  all  the  values  of  Y  that  belongs  to 
c  [-B-A,B-A]  such  that  Rl(Y)  is  between  two  given  levels 
c  LEVA  and  LEVS.  MERGESORT  is  called  to  obtain  the  new  lists 
c  of  zeros  and  teatpoints. 

ELSE 

LEVA=I - 1 
LEVB=I 


CALL  NERGEANOSORT  ( LEVA, LEVB, ZERO, TESTPT, 

*  MZERO , NTESTPT , NUMOF  ZERO ) 

NUM*NUMOF ZERO(  I  - 1  )+NUMOFZERO<  I ) 

DO  40  J«1,NUM 

IF  ( (NTESTPT ( J , 1 } . GT . LEV( I  - 1 ) ) .AMO . 

*  (MTESTPTC J, 1 ).LT.LEV( I )))  THEM 
IF  (J.EQ.1)  THEN 

ARS<L,I,K)"NZERO(J)*2 

ARSa,t,K+1)*ttZERO(J) 

urite<*, 100)  ars(l,i,k),ars(l,i,k«-1),k 

K=K*2 

ELSE 

ARS(L, I ,K)=MZERO( J- 1 ) 

ARS(L, I ,K+1 )=HZERO( J) 
write<*,100)  ars(l,i,k),ars(l,i,k+1),k 
K-K+2 
END  IF 
END  IF 

40  CONTINUE 

IF  ((MTESTPT(NUM,2) .GT.LEV( I  - 1 )) .AND. 

*  (MTESTPT(NUN,2).LT.LEV(I)))  THEN 
ARS(L, I ,K)=MZERO(NUM) 

ARS(L, I ,k+1 )=HZERO(NUM)-3 

write<*, 100)  ars( l, i , k),ars( l , i ,k+1),k 

K=K*2 

END  IF 
END  IF 
URITE(*,*) 

UR1TE(*,*) 

10  CONTINUE 

100  FORMAT(2E15 .5, 14) 

RETURN 

END 


c . copy  ars  into  ar . 

c  The  new  contents  stored  in  ARS  is  loaded  into  AR. 
c  If  TEST=0  then  after  each  iteration  the  net*  partition 
c  is  copied  into  the  file  "start2". 
c  If  TEST=1  or  2  then  after  each  iteration  the  new 
c  partition  is  copied  into  "start  1"  or  "start2" 

SUBROUTINE  UPDATE  (ARS.AR) 

INCLUDE  'define. f' 

INTEGER  DEL, I , J,K,R£F 

DOUBLE  PRECISION  ARS(NNS , NN , NMAX ) , AR (NNS , NN , NMAX ) 


IF  (NS.LE.2)  GO  TO  10 

URITE(*,*)  '  ERROR  NEED  TO  CHANGE  "UPDATE"  IN  ROOTFB.F. ' 

10  CONTINUE 
I =SNUMB 

DO  20  J=1,N 
DO  30  K=1,MAX 

ARC  I , J ,K)*ARS( I , J,K) 

30  CONTINUE 

20  CONTINUE 

IF  (CONO.EO.I)  THEN 

OPEN(UNfT>12, FILE"' Start 1 ' .STATUS" 'UNKNOWN ' ) 

REF*12 

IF  (LNAX)  THEN 

WR1 TE(*,*)  'THE  NEW  PARTITION  IS  LOADED  INTO  "startl".' 
END  IF 


% 


ELSE 


0PEN!UNIT=99, FILE*'start2' , STATUS- 'UNKNOWN' ) 

R£F=99 

WRITE!*,*)  'THE  NEW  PARTITION  IS  LOADED  INTO  "start2" 

END  IF 
DEL2- 1 

DO  15  1*1, NS 
OO  25  J«1,N 

DO  35  K*1 ,MAX,2 

IF  !AR!I, J,K).NE.AR!I  , J.K+1))  THEN 
WRITE!REF,120)  AR( I , J,K),AR( 1 , J.K+1 > 
WRITE(REF, 130)  DEL, DEL 

120  FORMAT!F15.9,3X,F15.9) 

130  FORMAT (216) 

ENO  IF 

110  FORMAT! F 16. 1 1 ) 

35  CONTINUE 

25  CONTINUE 

15  CONTINUE 

CLOSE  (UN1T=REF) 

RETURN 

END 


c  This  routine  Mill  clear  the  folloMing  matrices: 
c  ZERO, TESTPT , MZERO, MTESTPT , NUMOFZERO . 

SUBROUTINE  CLEAR  (ZERO, TESTPT , MZERO, MTESTPT , NUMOFZERO) 

INCLUDE  'define. f' 

INTEGER  I,J,K 

DOUBLE  PRECISION  ZERO! ZN.ZMAX), TESTPT (ZN.ZMAX, INTER) 
INTEGER  NUMOFZERO(ZN) 

DOUBLE  PRECISION  MZERO(MMAX) .MTESTPT (MMAX.M INTER) 

DO  10  1=1, ZN 

DO  20  J=1 , ZMAX 

ZERO! I , J )=O.ODO 

20  CONTINUE 

10  CONTINUE 

DO  100  1=1, ZN 

DO  200  J=1 ,ZMAX 

DO  300  K=1, INTER 

TESTPT! I , J,K)=0.000 
300  CONTINUE 

200  CONTINUE 

100  CONTINUE 

DO  15  1=1, ZN 

NUMOFZERO! I )=0. 000 
15  CONTINUE 

DO  11  1*1, MMAX 

DO  21  J=1 .MINTER 

MTESTPT! I, J)=0. 000 

21  CONTINUE 

11  CONTINUE 
DO  14  1=1 ,MMAX 

MZERO! I >=0.000 
14  CONTINUE 

RETURN 
ENO 


C . Ut(lNOX,AR) . 

c  Find  new  partition  and  store  it  in  the  matrix  AR. 
c  This  routine  will  be  called  in  the  MAIN  PROGRAM. 


SUBROUTINE  UKINDX.AR) 

INCLUDE  'define. f» 

INTEGER  INDX 

DOUBLE  PRECISION  AR(NNS,NN,NMAX),ARS(NNS,NN,NHAX) 
DOUBLE  PRECISION  ZERO(ZN,ZMAX),TESTPT(ZN,ZMAX, INTER) 
INTEGER  NUMOFZERO(ZN) 

DOUBLE  PRECISION  MZERO(MMAX)  ,MTESTPTO#tAX,M  INTER) 
EXTERNAL  C l earAR , CLEAR , L I STOFZERO, NEWPART , UPDATE 

CALL  Cl earAR  (ARS) 

CALL  CLEAR  ( ZERO, TESTPT,MZERQ,MTESTPT,NUMOF ZERO) 

CALL  LEVELS(LEV) 

CALL  LISTOFZERO  (ZERO, NUMOFZERO.TESTPT, LEV) 

CALL  NEWPART  (ZERO, TESTPT, I NDX.NUMOFZERO, ARS, LEV) 
CALL  UPDATE  (ARS,AR) 

RETURN 

END 


c 


JOHN.  F 


c  INITIALIZATION  OF  ARRAYS  USING  SUBROUTINE:  INITDI 
c 

c  ****  THIS  ROUTINE  IS  CHANGED  DEPENDING  ON  THE  EXAMPLE  CHOSEN 

c 

subroutine  initdi 
c 

integer  dist 
double  precision  a,b 
c 

c  Initialize  common  blocks  for  distributions, 
c 

include  'define. f' 
c 

c  Read  data  from  file 
c 

opentuni t=3, status* ' old' ,name=' distr.dat ' ) 
do  10  dist  =  azero,ns 

read  (3,*)  Uinfoldist, loner),  Uinfo<dist, upper) 
a=Uinfoldist, lower) 
b=Uinfoldist, upper) 

Winfo(dist,mean1 )=(b+a)/2.0d0 

Uinfoldist, ms)  =  Ub-a)**2)/12.0d0  ♦  Uinfoldist, meanl )**2 
10  continue 
close  (unit=3) 


c  Print  info,  and  compute  heights  for  uniform  distributions. 

c 

print  *,  'Underlying-distribution  information:' 
c 

do  20  dist  =  azero,ns 
Uinfoldist, height) 

&  =  1 .dO/lUinfoldist,upper)-Uinfoldist, tower)) 

IF  (OIST.EQ.0)  THEN 

UR  I TE(*, 1 1 )  Uinfoldist, lower), 

&  Uinfoldist, upper),  Uinfoldist, height) 

11  FORMAT l ’  X  : ' , F8.5, ’ , ' , F8.5, ' , ' , F8.5) 

ELSE 

URI TEl*, 12)  dist,  Uinfoldist, lower), 

&  Uinfoldist, upper),  Uinfoldist, height) 

12  FORMAT l'Z',I1,'  : ' , F8.5, ' , ' ,F8.5, ' , ' ,F8.5) 

END  IF 

20  cont i nue 

c 

c  print  *,  'Computed  observation-distribution  information:' 
c 

do  30  dist  =  1 ,ns 

Yinfoldist, lower)  =  Uinfoldist, lower)  ♦  Uinfolazero, lower) 
Yinfoldist, upper)  =  Uinfoldist, upper)  ♦  Uinfolazero, upper) 
c  print  *,  dist,  '  :  ',  Yinfoldist, lower),  Yinfoldist, upper) 

30  continue 

c 

return 

end 


c 

c 

c 

c 

c 

c 

c 


CINTGD  function 

This  function  returns  the  required  integrand  for  computing 

probl  Y_s1  in  lal.bl]  |  Y_s 2  =  yy  )  *  f_<Y_s2>(y) 

l set  igtype=cndprb*1) 


**** 


(set  igtype=cndexp=2) 


c  Et  x  I  Y_s2  =  yy  1  *  f_CY_s2)(yy) 
c 

c  density  for  r_s2  (set  igtype=convol=3) 

c 

double  precision  function  cintgd(x) 
include  'define. ♦' 

double  precision  densty,  PW,  x,  z,  l,  u 
c 

c  Integrand  for  computing  density  of  Y_s2,  which  is  the  convolution 
c  of  the  density  of  X  =  U_azero  with  the  density  of  W_s2. 
c 

l  =  yy-x 

z  =  densty(azero.x)  *  densty(s2,l) 
c 

goto  (  110,  120,  130  ),  igtype 
c 

c  Integrand  for  computing  prob(  Y_s1  in  (a1,b1]  |  Y_s2  =  yy  ) 
c 

110  l  =  lll-x 

u  =  ull-x 
z  =  z  *  PU(s1,l,u) 
goto  130 
c 

c  Integrand  for  computing  Et  X  |  Y_s2  =  yy  ] 
c 

120  z  =  z  *  x 

c 

130  cintgd  =  z 

c 

return 

end 


c  FUNCTION  UINTGO,  set  up  integrand  for  nonconditional  "stuff". 

c 

c  ****  THIS  ROUTINE  IS  CHANGED  DEPENDING  ON  THE  EXAMPLE  CHOSEN  *»** 

C 

c 

c  This  function  returns  the  required  integrand  for  computing 
c 

c  EC  X  I_( l 11 ,ul1] (Y_s1 )  ]  (set  igtype=corr=1) 

c 

c  prob(  Y_s1  in  (lll.ulll  ,  Y_s2  in  (U2,ul21  )  (set  igtype=  jpr=2) 
c 

c  prob(  Y_s1  in  (lll.ulll  )  (set  igtype=unipr=3) 

c 

double  precision  function  uintgd(x) 
include  'define. f' 

double  precision  densty,  PU,  x,  z, l,u, l l ,ul,a,b,zplus 
logical  novoid 
external  inters 
c 

c  Integrand  for  computing  prob(  Y_s1  in  (lll.ulll  ). 

c 

l  =  lll-x 
u  =  ull-x 

z  =  densty (azero.x)  *  PU(s1,l,u) 

c 

goto  (  110,  120,  130  ),  igtype 

c 

c  Integrand  for  computing  Et  X  I_(ll1,ul1J(Y_s1)  1 
c 

110  if  (Imax)  then 


so 


a=winfo(s1 , tower) 
b=winfo(s1,ipper) 
call  inters(a,b,l,u, It, ul, novoid) 
if  (novoid)  then 

zplus»  winfo(s1 .height)*. 5*(ut**2- l  l**2)*densty(azero,x) 

else 

zplus=0.0d0 
end  if 

else 

zplus  =  O.OdO 

end  i  f 

z  =  z  *  x  +  zplus 

goto  130 
c 

c  Integrand  fc r  computing  prob(  Y_s1  in  <111, uL 13  ,  Y_s2  in  (112, ul 2]  ) 
c 

120  l  =  l l 2 - x 

u  =  ul2-x 
i  -  z  *  PU(s2, l , u) 
c 

130  uintgd  =  z 

return 

end 


c  ****  THIS  ROUTINE  IS  CHANGEO  DEPENDING  ON  THE  EXAMPLE  CHOSEN  **** 

c 

c  Compute  densty(dist.x),  where  dist=0  implies  X 
c  and  dist>=1  implies  W_dist 
c 

double  precision  function  densty(dist.x) 
double  precision  x 
integer  dist 
include  'define.f' 
c 

densty  =  Uinfo(dist, height) 

return 

end 


c  ****  THIS  ROUTINE  IS  CHANGED  DEPENDING  ON  THE  EXAMPLE  CHOSEN  **** 

c 

c  Compute  prob(  W_dist  (a,bl  ) 

c 

double  precision  function  PU(dist,a,b) 

double  precision  a,  b,  l,  u,  llim,  ulim,  z 
integer  dist 

include  'define.f' 
c 

if  (  b  .It.  a  )  then 
z  =  0 
else 

llim  =  Winfo(dist, lower) 
ulim  =  Uinfo(di st, upper) 
if  (  b  .It.  llim  )  then 
z  a  0 

else  if  (  a  ,gt.  ulim  )  then 
z  =  0 

else 

if  (  a  .It.  llim  )  then 
l  -  llim 
else 
l  =  a 
endi  f 


£1 


if  (  b  .gt.  ulim  )  then 
g  =  ulim 
else 
u  =  b 
endif 

z  *  Winfo(dist,height)*(u  -  t) 

endif 
endif 
PW  =  z 
return 
end 


c 

c  Determine  the  intersection  of  [a,b]  and  [c,d] .  Set  novoid  =  .true, 
c  if  the  intersection  is  nonempty.  Set  novoid  =  .false,  otherwise, 
c  If  novoid  =  .true,  then  [l,u]  is  the  intersection  of  [a,b]  and  [c,d]. 
c 
c 

SUBROUTINE  inters(a,b,c,a, l ,u,novoid) 


double  precision  a,b,c,d,l,u 
logical  novoid 


c 


if  (b  .It.  a)  then 

novoid=  .false, 
else  if  (d  .It.  c)  then 
novoid=  .false, 
else  if  (  d  .It.  a  )  then 
novoid  =  .false, 
else  if  (  c  .gt.  b  )  then 
novoid  =  .false. 

else 


novo id 
if  (  d 

else 

endi  f 
if  (  c 

else 

endi  f 


=  .true. 

.It.  b  )  then 
u  =  d 

u  =  b 

-gt.  a  )  then 
l  =  c 

l  =  a 


endi  f 

return 

end 


c  This  routine  decodes  an  interval 

SUBROUTINE  FIXUP  (LlX,ULX,ll,UL,MAXlfMAXU) 

DOUBLE  PRECISION  CHECK 

DOUBLE  PRECISION  LLX,ULX,LL,UL,NAXU,HAXL 

CHECK=ULX-LLX 

c  Ca.+INF]  - >  la,a-3]  then  check=-3.0d0  NB:  ♦INF=NAXU 

IF  (CHECK. LT. -2. 600)  THEN 
UL=MAXU 
LL=LLX 

c  C-INF,b]  - >  [2+b,b]  then  check=-2.0d0  NB:  -iNF=MAXL 

ELSE  IF  (CHECK.LT. -1.600)  THEN 


UL=ULX 

LL=MAXL 

c  C- INF,+INF]  - >  10,-1]  then  checks- I.OdO 

ELSE  IF  (CHECK.LT. -.500)  THEN 
UL=MAXU 
LLsMAXL 

c  If  check  >0  then  we  have  a  finite  interval. 
ELSE 

UL=ULX 
LL=LLX 
END  IF 
RETURN 
END 


c 
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INITIALIZATION  OF  ARRAYS  USING  SUBROUTINE:  INITOI 
subroutine  initdi 
integer  dist 

double  precision  a,b, First_Mement,Sec_Moment,  Imse, alpha 
externa t  F i rs t_Mo»ent , Sec_Moment 

Initialize  common  blocks  for  distributions. 

include  'define. f' 

Read  data  from  file 

open(unit=3,status='old' ,name='distr.dat' ) 
do  10  dist  =  azero.ns 

read  (3,*)  WinfoCdist,  lower),  Winfo(dist, upper) 
a=Winfo(dist, lower) 
b=Winfo(di st, upper) 

Winfo(dist,mean1 )=First_Moment(b,dist)-Fi rst_Moment(a,dist) 
Uinfo(dist,ms)=Sec_Moment(b,dist)-Sec_Moment(a,dist)-Uinfo(dist,mean1)**2 
10  continue 
close  (unit=3) 


a=Winfo(azero,ms) 
b=Uinfo(1 ,ms) 

write(*,*)  a,'  and  ',b,'  are  a  and  b' 

alpha=a/(a+b) 

lmse=a*(1 .OdO-alpha) 

wri te(*,*) 

write(*,*)  ' LMSE  is  equal:  ',tmse 
wri te(*,*) 


c  Print  info,  and  compute  heights  for  uniform  distributions. 

c 

print  *,  'Underlying-distribution  information:' 
c 

do  20  dist  =  azero.ns 

WRITE(*,11)  dist,  Winfotdist, lower), 

&  Winfofdist, upper) 

11  FORMAT ( 14, ' : ' , F8.5, ' , ' , F8.5) 

20  continue 

c 

c  print  *,  'Computed  observation-distribution  information:' 
c 

do  30  dist  *  1,ns 

Yinfo(dist, lower)  =  Uinfo(dist, lower)  ♦  Winfo(azero, lower) 
YinfoCdi st, upper)  =  Winfofdist, upper)  +  Uinfo(azero, upper) 
URITEC*, 12)  dist,  Yinfofdist, lower), 

S  Yinfo(dist, upper) 

12  FORMAT ('Y',I4,':',F8.5,',',F8.5) 

30  continue 


return 

end 


c  CINTGO  function 


c  This  function  returns  the  required  integrand  for  computing 


c  prob(  Y_s1  in  (al.bll  |  Y_s2  =  yy  ) 

c 

c 

c  Et  X  I  Y_s2  *  yy  )  *  f_<Y_s2)<yy) 
c 

c  density  fop  Y_a 2 

c 

double  precision  function  cintgd(x) 
include  'define. f* 

double  precision  densty,  PW,  x,  z,  l,  u 


*  f_(Y_s2Xy) 

(set  igtype*cndprb=1) 

(set  igtype*cndexp=2) 

(set  igtypesconvoi=3) 


c  Integrand  for  computing  density  of  Y_s2;  which  is  the  convolution 
c  of  the  density  of  X  =  U_azero  with  the  density  of  U_s2. 
c 

l  =  yy-x 

z  *  densty(azero.x)  *  densty(s2,l) 
c 

goto  (  110,  120,  130  ),  igtype 
c 

c  Integrand  for  computing  prob(  Y_s1  in  (al.bl]  j  Y_s2  =  yy  > 
c 

110  l  =  lll-x 

u  =  ull-x 
z  *  z  *  PW(s1 , l ,u) 
goto  130 

c 

c  Integrand  for  computing  Et  X  |  Y_s2  =  yy  ] 

c 

120  z  =  z  *  x 


130  cintgd  =  z 

return 

end 


c  FUNCTIOM  UINTGD,  set  up  integrand  for  nonconditional  “stuff". 

c 

c  This  function  returns  the  required  integrand  for  computing 

c 

c  Et  X  I_< 111 ,ul1] (Y_s1 )  ]  (set  igtype=corr=1) 

c 

c  prob(  Y_s1  in  (lll.ull]  ,  Y_s2  in  (112, ul2]  )  (set  igtype=jpr=2) 
c 

c  prob(  Y_s1  in  (lll.ull]  )  (set  igtype=unipr=3) 

c 

double  precision  function  uintgd(x) 
include  'define. f' 

double  precision  densty,  PW,  x,  z, l,u, ll,ul,a,b,zplus,First_Moment 
logical  novoid 

external  inters,First_Moment 

c 

c  Integrand  for  computing  prob(  Y_s1  in  (lll.ull]  ). 
c 

l  *  lll-x 
u  =  ull-x 

z  =  densty(azero.x)  *  PW(s1,l,u) 
c 

goto  (  110,  120,  130  ),  igtype 
c 

c  Integrand  for  computing  Et  X  I_(ll1,ul1](Y_s1)  ] 


5S 


110  if  (tmax)  then 
a*winfo(s1 , lower) 
b*winfo(a1, upper) 
call  inters(a,b, l,u, ll.ul.novoid) 
if  (novoid)  then 

Zplu#*(First_Ho«ent(b,s1 )-First_Mo*ent(a,s1))*densty(azero,x) 

else 

zplus a0.0d0 
end  if 

else 

zplus  =  O.OdO 

end  i  f 

z  =  z  *  x  +  zplus 

goto  130 

Integrand  for  computing  prob<  Y_s1  in  (lll.ull)  ,  Y_s2  in  (112, ul2]  ) 

120  l  =  U2-x 

u  =  ul2-x 
z  =  z  *  PU(s2, l,u) 

130  uintgd  =  z 

return 

end 


Compute  densty(dist,x),  where  dist=0  implies  X 
and  dist>=1  implies  W_dist 

double  precision  function  densty(dist.x) 

include  'define. f' 

double  precision  x,X01aaf ,PI , alpha, b,v, heights 
integer  dist 
external  xOlaaf 

Pl-XOIaaf (x) 
a=Uinfo(dist, lower) 
b=Uinfo(di st, upper) 
v=(3.0d0*PI )/(2.0d0*b) 
alpha=(6.0dO*PI)/((15.0dO*PH-8.0dO)*b) 
heights2  1.d0/(b-a) 

If  (dist.eq.1)  then 

densty=  alpha*(5.0d0/4.0d0  -  cos(v*x)) 
densty  =  heights 
else  if  (dist.eq.azero)  then 
densty  *  heights 

densty=  alpha*(5.0d0/4.0d0  -  cos(v*x)) 
else 

densty  =  heights 

densty=  alpha*(5.0d0/4.0d0  -  cos(v*x)) 
end  if 

return 

end 

Compute  prob(  U_dist  in  (a,W  ) 

double  precision  function  PU(dist,a,b) 

double  precision  a,  b,  l,  u,  llim,  ulim,  z,Distr_Func 
integer  dist 

include  'define. f' 


external  Distr  Func 


c 


if  (  b  .It.  a  )  then 
z  *  0 
else 

Him  =  Winfofdist, lower) 
ulim  *  Uinfo(dist,<4>per) 
if  (  b  .It.  Him  )  then 
z  *  0 

else  if  <  a  .gt.  ulim  )  then 
z  =  0 

else 

if  (  a  .It.  Him  )  then 
l  =  llim 
else 


l  =  a 
endif 

if  (  b  .gt.  ulim  )  then 
u  =  ulim 
else 
u  =  b 
endif 


z=  Distr_Func(u,dist)-Distr_Func( l ,dist) 

endi  f 
endif 
PW  =  z 
return 
end 


c 

c  Determine  the  intersection  of  [a,b]  and  Cc,d].  Set  novoid  =  .true, 
c  if  the  intersection  is  nonempty.  Set  novoid  *  .false,  otherwise, 
c  If  novoid  =  .true,  then  [l,u]  is  the  intersection  of  [a,b]  and  [c,d] . 

c 

c 

SUBROUTINE  inters<a,b,c,d, l ,u, novoid) 


double  precision  a,b,c,d,l,u 
logical  novoid 


c 


if  (b  .It.  a)  then 

novoid=  .false, 
else  if  (d  .It.  c)  then 
novoid=  .false, 
else  if  (  d  .It.  a  )  then 
novoid  =  .false, 
else  if  (  c  .gt.  b  )  then 
novoid  =  .false. 


else 


novoid  =  .true, 
if  (  d  .It.  b  ) 
u  =  d 

else 

u  *  b 


endif 

if  (  c  .gt.  a  ) 
l  *  c 

else 

l  =  a 


endif 


then 


then 


endif 


sn 


return 

end 


SUBROUTINE  FIXUP  (LLX,ULX,LL,UL,MAXL,MAXU) 

DOUBLE  PRECISION  CHECK 

DOUBLE  PRECISION  LLX,ULX,LL,UL,MAXU,HAXL 

CHECK=ULX-LLX 

c  [a,+lNF]  - >  Ca,a-3]  then  check=-3.0d0  NB:  +INF=KAXU 

IF  (CHECK. LT . -2.600)  THEN 
UL=MAXU 
LL=LLX 

c  [■ INF ,b]  - >  C2+b,b]  then  check=-2.0d0  NB:  -INF=MAXL 

ELSE  IF  (CHECK.LT. -1.600)  THEN 
UL=ULX 

LL=MAXL 

c  [-INF.+INF]  - >  [0,-1]  then  check=-1.0d0 

ELSE  IF  (CHECK.LT. -.500)  THEN 
UL=MAXU 
LL=MAXL 

c  If  check  >0  then  we  have  a  finite  interval. 

ELSE 

UL=ULX 
LL=LLX 
END  IF 
RETURN 
END 


double  precision  function  Sec_Moment  (X, point) 

include  'define. V 
integer  point 

double  precision  x, alpha, v,b, PI ,x01aaf,a,heights 
external  xOlaaf 

pi=X01aaf(x) 
a=Uinfo(point, lower) 
b=Winfo(point, upper) 
v=(3*PI )/(2*b) 

alpha=(6.0d0*PI )/((15.0d0*PI+-8.0d0)*b) 
heights®  1.d0/(b-a) 

If  (point. eq.azero)  then 
c  Sec_Moment=((X**3)*heights)/3.0d0 

Sec_Moment=( ( 5 . OdO/ 1 2 . OdO )*X**3- ( 1 /v)**3*( v*x*cos( v*x ) -  2 . OdO*s i n<  v*x ) 
&  + v**2*x**2 . OdO*a i n( v*x ) ) )*A  t  pha 

else  if  (point. eq.1)  then 
c  Sec_Moment=((X**3)*heights)/3.0d0 

Sec_Moment=((5.0d0/12.0d0)*X**3-(1/v)**3*(v*x*cos(v*x)-2.0d0*sin(v*x) 
&  "  ♦v**2*x**2.0d0*sin(v*x)))*Atpha 

else 

c  Sec  Noment»((X**3)*heights)/3.0d0 

Sec~Haeent*((5.0d0/12.0d0)*X**3-(Vv)**3*(v*x*cos(v*x)-2.0d0*sin(v*x> 
&  “  +v**2*x**2.0d0*sin(v*x)))*Alpha 

end  if 

return 

end 
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double  precision  function  Firstjtaaent  (X, point) 

include  'define. f' 
integer  point 

doti>le  precision  x, alpha, v,b,PI, xOlaaf , a, heights 
external  xOlaaf 

pi*X01aaf (x) 
a=uinfo(point, lower) 
b=Uinfo( point, upper) 
v*(3.0dO*P! )/(2.0d0*b) 
alpha=(6.0dO*PI )/((15.0dO*PI+8.0d0)*b) 
heights*  1.d0/(b-a) 

If  (point. eq.atero)  then 
First  Moment=(X**2)*.50d0*heights 

First~Moment=((5.0d0/8.0d0)*x“2  -(1/v)“2*(cos(v*x)+v*x*sin(v*x))) 
&  "  ‘Alpha 

else  if  (point. eq.1)  then 

First  Moment=(X“2)*.50d0*heights 

First“Homent=((5.0d0/8.0d0)*x“2  -(1/v)**2*(cos(v*x)+v*x*sin(v*x))) 

&  ‘Alpha 

else 

First  Homent=(X“2)*.S0d0*heights 

F i  rst~Moment=( (5 . OdO/8. OdO )*x“2  -(1/v)“2*(cos(v*x)+v‘x*sin(v*x))) 
&  ‘Alpha 

end  i  f 

return 

end 


double  precision  function  Distr_Func  (X, point) 

include  'define. f' 
integer  point 

double  precision  x, alpha, v.b.PI , xOlaaf , a, he:  ghts 
external  xOlaaf 

pi*X01aaf(x) 
a=Winfo( point, lower) 
b*Uinfo( point, upper) 
v=(3.0dO*PI )/(2.0d0*b) 
alpha=(6.0dO*PI )/((15.0d0*PI+8.0d0)*b) 
heights*  1.d0/(b-a) 

If  (point. eq.azero)  then 
Oistr_Func=(X)*heights 

Distr_Func*((5.0d0/4.0d0)*x  -(1/v)*sin(v*x))*Alpha 
else  if  (point. eq.1)  then 
Oistr_Func=((5.OdO/4.0d0)*x  -(1/v)‘sin(v*x))*Alpha 
Oistr_Func=(X)*heights 
else 

Distr_Func*(X)*heights 

0istr_Func*((5.0d0/4.0d0)‘x  -<1/v)*sin<v*x))*Alpha 
end  if 


return 

end 


c 


CONO.  F 


DOUBLE  PRECISION  FUNCTION  CONV  (IHD,Y) 

INCLUDE  'define. f' 

INCLUDE  'nefl.f' 

INTEGER  IND 

DOUBLE  PRECISION  Y,L,U,A,B,C,D,RES, CINTGO 
LOGICAL  NOVOID 
EXTERNAL  CINTGO, INTERS 

A=WINFO(AZERO, LOWER) 

8=WINF0(AZERO, UPPER) 

S2=IND 

YY=Y 


IF  (YY.GT .YINFO(S2, LOWER)  .AND.  YY.LT .YINFO(S2, UPPER))  THEN 
C=YY-WINFO(S2, UPPER) 

D=YY-WINF0(S2, LOWER) 

CALL  INTERS(A,B,C,0,L,U, NOVOID) 

IF  (NOVOID)  THEN 
IGTYPE=CONVOL 


EPSABS=  Id- 10 
EPSREL=1d-8 
I FAIL=- 1 

CALL  D01AJF(CINTGD,L,U,EPSABS,EPSREL,RES 
*  ,ABSERR,W,LW, IW,LIW,IFAIL) 

IF  (IFAIL.NE.O)  THEN 

WRITE(*,*)  'IN  CONV' 

WRITE!*,*)  'LOWER  LIMIT:  ',L 
WRITE(*,*)  'UPPER  LIMIT:  ',U 

END  IF 

ELSE 

RES=0.000 
END  IF 

ELSE 

RES*0.000 


END  IF 
CONV=RES 
RETURN 
END 


OOUBLE  PRECISION  FUNCTION  EC(CONDI.Y) 

INCLUDE  'define. f' 

INCLUDE  ' nag . f ' 

INTEGER  CONDI 

DOUBLE  PRECISION  Y,L,U,A,B,C,D,RES,CONV,C!NTGD 

LOGICAL  NOVOID 

EXTERNAL  CINTGO, CONV, INTERS 

A=WINFO(AZERO, LOWER) 

B=WINFO(AZERO, UPPER) 

S2=CONDI 

YY=Y 

IF  (YY.GT. Y1NF0(S2, LOWER)  .AND.  YY.LT. YINFO(S2, UPPER))  THEN 
IF  (LMAX)  THEN 
RES-Y 

ELSE 

C»YY-WINFO(S2, UPPER) 

D=YY-WINFO(S2, LOWER) 

CALL  INTERS(A,B,C,D,L,U, NOVOID) 

IF  (NOVOID)  THEN 


(oO 


I GTYPE=CNOEXP 


* 


ELSE 

END  IF 
EC=RES 
RETURN 
END 


EPSABS*  Id- 10 
EPSR£L*1d-8 
I FAIL--1 

CALL  D01AJF(CINTGD,L,U, EPSABS, EPSREL, RES 
,ABSERR,U,LW, IW, LIU, I  FAIL) 

IF  (IFA1L.HE.0)  THEN 

WRITE(*,*)  'IN  EC.' 

WRITEf*,*)  'LOWER  LIMIT:  ',L 
WRITE<*,*)  'UPPER  LIMIT:  ',U 

END  IF 

RES=R£S/C0NV(S2,YY) 

ELSE 

RES=0.000 
END  IF 
END  IF 

RES-O.ODO 


DOUBLE  PRECISION  FUNCTION  FUN4  (LLX,ULX,Y, INO, CONDI ) 

INCLUDE  'define. f' 

INCLUDE  ' nag . f ' 

INTEGER  I ND, COND i 

DOUBLE  PRECISION  Y,L,U,A,B,C,0,RES,CONV,CINTGD 
DOUBLE  PRECISION  LLX,ULX,MAXL,MAXU 
LOGICAL  NOVO ID 

EXTERNAL  CINTGD, CONV , INTERS, FIXUP 


A*WINFO(AZERO, LOWER) 

B=UINFO(AZERO,UPPER) 

SI -I  NO 

SZ=CONOI 

YY=Y 

MAXU=YINF0(S2, UPPER) 

MAXL=YINF0(S2, LOWER) 

CALL  FIXUP  (LLX,ULX,LL1,UL1 ,MAXL,MAXU) 

IF  (YY.GT.YINF0(S2, LOWER)  .ANO.  YY.LT.YINF0<S2, UPPER))  THEN 
C=YY-UINF0(S2, UPPER) 

D=YY-WINF0(S2, LOWER) 

CALL  INTERS  (A,B,C,D,L,U.NOVOID) 

IF  (NOVOID)  THEN 
IGTYPE=CNDPRB 

EPSABS=  Id- 10 
EPSREL=1d-8 
EPSREL=1d-6 
I  FA  I L  =  - 1 

CALL  D01AJF<CINTGD,L,U, EPSABS, EPSREL, RES 
*  ,ABSERR,U,LU, IW,LIW, IFAIL) 

IF  (IFAIL.NE.O)  THEN 

WRITE(*,*)  'IN  FUN4. ' 

WRITE<*,*)  'LOWER  LIMIT:  ',L 
WRITE(*,*)  'UPPER  LIMIT:  ',U 
URITEC*,*)  'THE  INO  and  COND  ARE:  '.IND.CONO 
END  IF 

RES=RES/C0NV(S2,YY) 

ELSE 


ELSE 


RES»0.000 
END  IF 

R£S*0.000 

END  IF 
FUN4=RES 
RETURN 
END 


c 


UNCONO . F 


DOUBLE  PRECISION  FUNCTION  PG  (LLX.ULX, INDEX) 

INCLUDE  'rwg.f ' 

INCLUDE  'defi ne.f' 

INTEGER  INDEX 

OOUBLE  PRECISION  LLX,ULX,L,U,R£S,MAXL,MAXU 
DOUBLE  PRECISION  UINTGD, A, B,C,D 
EXTERNAL  UINTGD, FIXUP, INTERS 
LOGICAL  NOVO 10 

IF  (LMAX)  THEN 
S1=COND 
ELSE 

SI = INDEX 
END  IF 

MAXU=YINFO(S1,UPPER) 

MAXL=YINFO(S1 , LOWER) 

CALL  FIXUP  (LLX.ULX, LL1 ,UL1 .MAXL.MAXU) 
A=WINFO(AZERO, LOWER) 

B=UINFO(AZERO, UPPER) 

C=LLl  -  WINFO(S1, UPPER) 

D=UL1  -  WINFO(S1, LOWER) 

CALL  INTERS(A,B,C,D,L,U, NOVOID) 

IGTYPE=UNIPR 
I FAI L=- 1 
EPSABS=1D- 10 
EPSREL=1 .00-08 


IF  (NOVOID)  THEN 

CALL  D01AJF(UINTGD,L,U,EPSABS,EPSREL,RES,ABSERR, 
*  W,LW, IW, LIW, I  FAIL) 

IF  ( I FAIL.NE.O)  THEN 


WRITE(*,*) 

'IN  PG' 

WRITE(*,*) 

'LOWER  LIMIT:  ',L 

WRITE(*,*) 

'UPPER  LIMIT:  ',U 

END  IF 
ELSE 

RES=0.000 

END  IF 
PG=RES 
RETURN 
END 

DOUBLE  PRECISION  FUNCTION  JPG(LLX1 ,ULX1 ,LLX2,ULX2, IND1 , IND2) 

INCLUOE  ' nag . f ' 

INCLUOE  'define. f' 

INTEGER  I ND 1 , I ND2 

DOUBLE  PRECISION  LLX2.ULX2.LLX1 ,ULX1 
DOUBLE  PRECISION  UINTGD, L.U, RES, HAXL,NAXU 
DOUBLE  PRECISION  A,B,C,D, TENPt ,TEHPH 
EXTERNAL  UINTGD, FIXUP, INTERS 
LOGICAL  NOVOID 

S1=IND1 

S2-IND2 

A=WINFO(AZERO, LOWER) 

B=WINFO(AZERO, UPPER) 

NAXL=YINFO(S1 .LOWER) 

MAXU=YINFO(S1 .UPPER) 

CALL  FIXUP  (LLX1 ,ULX1 ,LL1 , UL1 ,HAXL,MAXU) 

C=  LL1  -  WINFO(S1, UPPER) 


0*  UL1  -  WINFO(S1, LOWER) 

CALL  INT£RS(A,B,C,D,TEMPL,TENPH, NOVOID) 

IF  (MOVOID)  THEN 
MAXL=YINF0(S2, LOWER) 

MAXU*YINF0(S2, UPPER) 

CALL  FIXUP  (LLX2.ULX2, LL2,UL2,MAXL,MAXU) 

C*  LL2  -  WINF0<S2, UPPER) 

0^  UL2  -  WINF0(S2, LOWER) 

CALL  IMTERS(TEMPL,TEMPH,C,D,L,U, NOVOID) 
c  write(*,*)  L,U 
IF  (NOVOID)  THEN 
IGTYPE=JPR 
I  FAIL*- 1 
EPSABS=1D-10 
EPSREL=1. 00-06 
c  EPSR£L=1 .00-08 

CALL  D01AJF(UINTGD,L(U,EPSABS,EPSREL,RES,ABSERR, 
*  W,LU,1U, LIU, [FAIL) 

IF  ( I FAIL.NE.O)  THEN 
WRITE!*,*)  'IN  JPG' 

WRITE(*,*)  'LOWER  LIMIT:  ',L 
WRITE(*,*)  'UPPER  LIMIT:  ',U 
END  IF 
JPG=RES 
ELSE 

JPG=0.000 
GOTO  10 
END  IF 
ELSE 

JPG=0.000 
END  IF 

10  CONTINUE 

RETURN 
END 


DOUBLE  PRECISION  FUNCTION  MEANX  (LLX.ULX, INDEX) 

INCLUDE  ' nag . f ' 

INCLUDE  'define. f' 

INTEGER  INDEX 

DOUBLE  PRECISION  LLX,ULX,L,U,MAXL,MAXU,RES 
DOUBLE  PRECISION  UINTG0,A,B,C,0 
LOGICAL  NOVO ID 

EXTERNAL  UINTGO, FIXUP. INTERS 


IF  (LMAX)  THEN 
S1=COND 
ELSE 

S1=INDEX 
END  IF 

MAXU=YINFO(S1 .UPPER) 

MAXL=YINFO(S1, LOWER) 

CALL  FIXUP  (LLX.ULX, LL1,UL1,MAXL,MAXU) 
A=WlNFO(AZERO. LOWER) 

B=W I N  FO( AZERO , UPPER ) 

C=LL1  -  UINFO(S1, UPPER) 

0=UL1  -  WINFO(ST, LOWER) 
c  urite(*,*)  C,d  '  Is  c  and  d' 

CALL  IMTERS(A,B,C,D,L,U, NOVOID) 
c  urite(*,»)  L,U 
IGTYPE=CORR 
IFAIL=-1 


EPSABS-1D-10 
EPSREL=1. 00-08 


IF  (NOVOID)  THEN 

CALL  DOIAJFOJINTGD, L,U,EPSABS,EPSREL,r>ES,ABSERR, 
*  U,LU, IV,LIW, IFAIL) 

IF  (IFAIL. NE.O)  THEN 

WRITE!*, *>  'IN  NEANX' 

WR1TE(*,*)  'LOWER  LIMIT:  ',L 
WR!TE(*,*)  'UPPER  LIMIT:  '.U 

END  IF 
ELSE 

RES=0.000 

END  IF 
MEANX=RES 
RETURN 
END 
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On  the  Computation  of  Shot-Noise  Probability  Distributions 
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Abstract  —  A  general  method  that  does  not  use  numerical 
integration  is  presented  for  approximating  a  continuous  cumu¬ 
lative  distribution  function  (cdf)  from  its  characteristic  func¬ 
tion.  If  the  cdf  has  a  density  that  is  piecewise  smooth  except 
for  jump  discontinuities,  and  if  the  discontinuity  locations  can 
be  identified,  then  excellent  density  approximation  can  be  ob¬ 
tained  by  fitting  a  cubic  spline  to  the  cdf  approximation  and 
then  differentiating  the  spline.  Since  shot-noise  densities  typ¬ 
ically  contain  essential  discontinuities  as  well  as  impulses,  the 
method  is  adapted  to  handle  these  complications. 

Index  Terms  —  Filtered  point  process,  Poisson  process, 
Fourier  transform,  Fourier  series,  fast  Fourier  transform,  Gibbs 
phenomenon,  sampling  theorem,  spline. 

I.  Introduction 

Shot-noise  processes,  also  known  as  filtered  point  pro¬ 
cesses,  constitute  an  important  class  of  mat  hematical  mod¬ 
els  used  to  understand  physical  phenomena  ranging  from 
the  measurement  of  nerve  impulses  in  the  brain,  to  the 
formation  of  images  on  film  exposed  under  low-level  illu¬ 
mination,  to  the  electric  current  generated  by  photodiodes 
used  in  optical  communication  systems.  Hence,  it  is  un¬ 
fortunate  that  in  most  cases,  shot-noise  densities  must  be 
obtained  by  numerical  contour  integration  of  their  charac¬ 
teristic  functions  or  other  sophisticated  techniques. 

We  first  present  a  general  method  for  recovering  a  con¬ 
tinuous  cumulative  distribution  function  (cdf)  from  its 
characteristic  function  without  numerical  integration.  If 
the  cdf  has  a  density  /(y)  that  is  piecewise  smooth  except 
for  jump  discontinuities,  and  if  the  discontinuity  locations 
can  be  identified,  we  can  easily  recover  the  density  as  well 
by  fitting  a  cubic  spline  to  the  cdf  approximation  and  then 
differentiating  the  spline.  Unfortunately,  shot-noise  densi¬ 
ties  typically  not  only  have  an  impulse  at  the  origin  but 
also  satisfy  limy_o+/(y)  =  oo.  We  must  therefore  study 
the  characteristic  function  more  carefully  and  modify  the 
initial  method  to  account  for  these  complications. 

The  paper  is  organized  as  follows.  In  Section  II  we 
summarize  the  general  method  for  the  case  of  a  piecewise- 
smooth  density  with  jump  discontinuities.  In  Section  III 
we  introduce  our  shot-noise  model  and  summarize  prior 
work  on  recovering  shot-noise  densities.  In  Section  IV  we 
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analyze  shot-noise  characteristic  functions  and  moment¬ 
generating  functions.  In  Section  V  we  modify  the  method 
of  Section  II  to  handle  shot-noise  random  variables.  Sev¬ 
eral  examples  are  considered  and  density  plots  are  pre¬ 
sented.  Since  shot-noise  characteristic  functions  them¬ 
selves  can  also  be  difficult  to  compute,  we  consider  an  ap¬ 
proximation  scheme  in  Section  VI. 

II.  The  General  Method 

We  propose  the  following  general  approach  for  recover¬ 
ing  a  continuous  cdf  from  its  characteristic  function.  The 
approach  is  based  on  the  following  lemma,  which  is  proved 
in  the  Appendix.  We  then  discuss  a  method  for  recovering 
the  corresponding  density,  assuming  that  it  is  piecewise 
smooth. 

Lemma  l:  Let  g  be  a  finite  nonatomic  measure  on  IR 
with  characteristic  function 

d>(u;)  :=  f  e^uizdfi(x). 

J  -  OO 

Then 

OO 

/i(0,oo)  =  lim  Y]  bn^(mr/L),  (1) 

L— oo  * — ' 
r»  =  —  oo 

where 

(  1/2,  n  =  0, 

6„  =  <  —  j/nir,  n  =  odd,  (2) 

(  0,  n  =  even  ^  0. 

Furthermore,  for  0  <  L  <  oo,  the  error, 

OO 

bn<t>{nir/L)  -  /r(0  ,  oo), 

n=-oo 

is  upper  bounded  by  fx(—  oo,  —  L]  and  lower  bounded  by 
-H(L,o o). 

If  we  now  fix  any  y  G  IR  and  consider  the  measure 
fiy(C)  =  ji(C  +  y),  it  is  an  easy  corollary  that 

OO 

/i(y,  oo)  =  lim  Y'  bn$(mr / L)e~inTy^ L ,  (3) 

L  —  oo  * 

n  =  — oo 

with  the  corresponding  upper  and  lower  error  bounds  being 
n(—oc ,—L  +  y]  and  —  n(L  +  y,  oo).  If  we  restrict  attention 
to  y  in  a  subinterval  [Li ,  L-i)  C  (  —  L,  L),  then 

ji(— co.— L  +  y]  <  jr(— oo, -L  +  L2] 


(4) 


and 


-n(L  +  y,o c)  > -n(L  ±  L\,x).  (5) 

Note  that  if  /i(-x,0)  =  0,  then  the  upper  bound  is  zero, 
and  if  p(0.  x )  =  0,  then  the  lower  bound  is  zero.  To  obtain 
more  easily  computable  error  bounds,  we  assume  that 

M(s)  :  =  J  e3*  d^(x)  <  x,  for  all  s  €  IR 

Then  for  y  €  [/-i.  Io]  and  s  >  0,  we  have 

/i(L  +  y,  x)  <  n(L  +  L i,oc) 

<  f ^  l  ) 

where  /+(x)  :=  supJ>0[s.r  —  /\‘(s)]  and  I\(s)  :  =  In.W(s). 
Similarly. 

/i(-x,  -L  +  y]  <  /i(-x, -L  +  L->] 

<  (6) 

where  /_ ( j- )  :=  supJ>0[— sx  —  A'(—  s)].  It  is  well  known 
from  the  theory  of  large  deviations  that  I±(x)  is  convex, 
continuous,  and  nonnegative.  Furthermore,  the  suprema 
are  achieved  for  some  ,s  >  0  [5,  pp.  7-9]. 

Remark:  Since  <J>  is  the  characteristic  function  of  /i, 
equation  (3)  can  be  regarded  as  a  sampling  theorem  for 
the  approximately  bandlimited signal  $(-;).  If<J>  were  truly 
bandlimited  to  (/.t,  Lo];  i.e.,  if  /<(  — x,  Z-i]  =  fi(L 2.  x)  =  0. 
then  (3)  would  hold  exactly  for  finite  L  >  Li  -  L 1,  and 
2rr/(  Ln  —  L 1 )  would  be  the  Nyquist  sampling  period.  When 
is  not  truly  bandlimited,  taking  L  <  oc  will  introduce 
aliasing,  the  effect  of  which  can  be  bounded  as  discussed 
above. 

To  use  the  formula  in  (3).  set 


Theorem  (Birkhoff  and  de  Boor  [2]/-  Let  F  be  a  func¬ 
tion  that  is  four  times  continuously  differentiable  on  an 
interval  [a,  6].  If  Fm  is  a  sequence  of  interpolating  cubic 
splines  satisfying  the  additional  constraints  F^(a)  =  F'(a) 
and  C(6)  =  F'(6),  then  Fm.  F'm,  F".  and  F"'  converge 
uniformly  to  F.  F',  F".  and  F'",  respectively,  assuming 
that  as  m  —  x,  (t)  the  number  of  knots  of  Fm  tends  to 
infinity,  (11)  the  maximum  knot  spacing  tends  to  zero,  and 
(Hi)  the  ratios  of  the  knot  spacings  of  Fm  are  bounded. 

If  the  measure  n  has  a  piecewise-smooth  density,  we  ap¬ 
proximate  it  by  differentiating  (the  negative  of)  the  cubic 
spline  approximation  of  F[  v  between  the  knots.  To  ob¬ 
tain  a  good  approximation,  it  is  necessary  to  place  triple 
knots  at  the  density  jump  points,  double  knots  at  corners, 
and  single  knots  at  corners  in  the  derivative  of  the  density. 
This  will  become  clearer  in  the  examples. 


Example  l:  Let  p  be  the  probability  measure  whose  den¬ 
sity  f(y)  is  one-half  the  sum  of  the  uniform  density  on 
[—  5,  .5]  and  the  standard  normal  density.  The  correspond¬ 
ing  characteristic  function  is 


$(•*>)  =  2 


sin(-^/2) 


w/2 


Since  <$(*:)  decays  like  1/w,  and  since  6n  decays  like  1/n. 
we  see  that  for  fixed  L,  the  series  in  (7)  converges  uniformly 
to  a  continuous  function  as  N  — *  x,  and  thus  the  finite 
series  will  not  exhibit  Gibbs  phenomenon.  Of  course,  if  we 
were  to  differentiate  (7),  the  new  series  coefficients  would 
decay  like  1/n,  and  Gibbs  phenomenon  would  appear.  We 
now  take  L  =  6,  Li  =  0,  and  L2  =  3.  For  the  error  bounds, 
it  suffices  to  consider  (6)  since  the  density  is  symmetric. 
Now, 


M(s) 


n.s(y)  =  Y.  ^(rxir/L)e-^'L.  (7) 

n  =  -,V 

Now  let  m  :  =  0  and  cn  :  =  bn<P(mr/L)  for  n  ^  0.  Observe 
that  c_„  is  the  complex  conjugate  of  cn  and  that  c\  =  0 
for  V  even  We  can  therefore  write,  for  even  1 V, 

,v_  1 

FIM»)  =  ba*{Q)  +  2ReYcne-3nwy/L- 

n  =0 

If  we  let  y  =  k  Ay,  where  Ay  =  2L/A',  then  the  sampies 

F'l  s(k  Ay)  for  k  =  0 . V  —  1  can  be  computed  with  an 

.V-point  fast  Fourier  transform  (FFT).  We  then  propose 
to  fit  a  cubic  spline  [3]  to  the  samples  of  F[  v(/fc  Ay)  that 
lie  in  the  range  /,->].  This  approximation  is  motivated 
by  the  following  result 


It  is  easily  verified  numerically  that  the  maximum  value 
of  — s(— 3)  —  In  A/(  —  s)  is  greater  than  5.17  and  occurs  at 
s  =»  3.  The  bound  is  then  e~5  17  <  .006.  We  next  took 
.V  =  256  and  obtained  66  samples  of  F£  N  in  the  range 
[0,3.02].  We  first  used  30  uniformly  spaced  knots  in  the 
proper  subinterval  [  0001,3]  to  generate  a  cubic  spline  ap¬ 
proximation  of  F£  v  using  the  NAG  library  subroutines 
E02BAF  and  E02BCF.  In  Fig.  1  both  f(y)  and  minus 
the  derivative  of  the  spline  are  plotted.  The  oscillations 
at  1/2  are  not  due  to  the  Gibbs  phenomenon  associated 
with  Fourier  series.  Rather,  they  are  due  to  the  fact  that 
a  cubic  spline  with  distinct  knots  is  twice  continuously  dif¬ 
ferentiable,  while  limAf— 00  F£  N(y)  is  not  even  continuous 
at  y  =  1/2.  Based  on  Fig.  1,  we  fitted  a  new  spline  with 
only  8  knots,  including  a  triple  knot  at  1/2  to  handle  the 
discontinuity  in  the  density.  The  knot  sequence  was  .0001, 
5,  5,  .5.  1,  1.5,  2,  3.  The  resulting  approximation  to  /(y) 


is  shown  if  Fig.  2.  It  is  interesting  to  compare  Fig.  2  with 
the  straightforward  FFT-approximation  obtained  from 

!{y)  =  7T  f  Y’(w)e“-’u'vdu> 

™  J -oo 

JV-l 

as  —  y  g>(n  Ay)f-,Bd"»iw 
27T 

n  =  -(N-l) 

by  taking  Aj  =  jt/L  and  y  =  A:  Ay,  where  Ay  = 
(27r/.V)/Aw  =  2L/N .  With  these  substitutions, 

/(F  Ay)  st  =■  LiO)+2Re  V'dne--'^*"  ,  (8) 

where  do  =  0  and  d„  =  ^(mr/ L).  Using  L  =  6  and  N  = 
256.  we  obtained  the  graph  shown  in  Fig.  3,  which  clearly 
shows  the  Gibbs  phenomenon  associated  with  the  Fourier 
series  of  a  function  with  a  jump  discontinuity. 

III.  Shot-Noise  Models 

Consider  the  real-valued  process  {Zt}  given  by 

Z,  =  J^W.TVj+V',.  (9) 

V 

where  the  {Tu}  are  points  of  a  Poisson  process  with  non¬ 
negative  intensity  A( •),  {/!„}  is  an  independent,  identically 
distributed,  nonnegative  "gain”  sequence,  and  { V't }  is  a 
zero-mean  Gaussian  process.  We  assume  {T„},  and 

{ l j }  are  mutually  statistically  independent.  The  deter¬ 
ministic  function  h  is  the  system  impulse  response,  or  point 
spread  function,  depending  on  the  application. 

.4.  Applications  of  the  Model 

The  model  described  by  (9)  is  quite  general.  If  t  repre¬ 
sents  time,  {Z  i }  can  be  used  as  a  model  of  the  photoelectric 
current  generated  in  the  receiver  of  an  optical  communica¬ 
tion  system.  If  t  represents  two-dimensional  spatial  posi¬ 
tion.  then  {Z(}  can  be  used  as  a  model  of  image  formation 
on  film. 

/  Optical  Communication  Systems 

In  optical  communication  systems,  Zt  models  the  elec¬ 
tric  current  delivered  by  an  amplifier  whose  input  signal  is 
the  output  of  an  avalanche  photodiode.  In  this  situation 
the  {Tj,}  denote  the  times  at  which  arriving  photons  are 
detected  in  the  receiver,  the  {/!„}  are  the  avalanche  gains, 
h  is  the  impulse  response  of  the  amplifier,  and  {Vj}  mod¬ 
els  thermal  noise  in  the  amplifier.  The  intensity  A(  )  of  the 
Poisson  point  process  is  determined  by  the  brightness  of 
the  fight  falling  on  the  photodiode. 

For  high  light  levels,  Z(  is  approximately  normal  (even 
if  no  Gaussian  noise  is  added)  (18,  19].  For  low  light  levels, 
the  Gaussian  approximation  would  not  be  so  accurate. 


2.  Image  Detection 

Consider  x-ray  images  on  film  with  point  spread  function 
/i(t ,  r ) ,  1 ,  t-  G  IR2.  Photons  striking  the  film  cause  “smears” 
instead  of  well-defined  points.  The  response  of  the  film  at 
a  point  t  €  IR2  due  to  a  photon  arriving  at  a  point  Tv  is 
Avh(t.Ti,).  This  gives  rise  to  the  model 

zt  =  r*,), 

u 

which  is  a  special  case  of  (9)  in  which  the  Gaussian  noise 
is  absent. 

B.  Summary  of  Prior  Work 

The  most  basic  shot-noise  statistics,  namely  the  mean 
and  variance,  were  reported  by  Campbell  in  1909  [6,  7]. 
Shot  noise  was  also  investigated  by  Schottky  in  his  1918 
paper  on  spontaneous  current  fluctuations  in  electric  con¬ 
ductors  [21].  In  1944-45,  Rice  [19]  gave  an  extensive  anal¬ 
ysis  of  shot  noise  when  the  underlying  Poisson  process  has 
a  constant  intensity.  In  particular,  he  showed  that  as  the 
intensity  tends  to  infinity,  the  probability  distribution  of 
the  shot  noise  tends  to  a  normal  distribution.  In  1971  Pa- 
poulis  [18],  considering  underlying  Poisson  processes  with 
time-varying  intensity,  gave  numerical  bounds  on  the  dif¬ 
ference  between  the  true  shot-noise  distribution  and  the 
Gaussian  approximation.  The  first  nonasymptotic  results 
concerning  the  cumulative  distribution  of  shot  noise  ap¬ 
peared  in  the  1960  paper  by  Gilbert  and  Poliak  [10].  For 
an  underlying  Poisson  process  with  constant  intensity,  they 
derived  an  integral  equation  satisfied  by  the  shot-noise  dis¬ 
tribution.  They  were  able  to  solve  this  integral  equation 
for  some  special  cases  of  the  shot- noise  impulse  response. 
Progress  on  the  numerical  computation  of  the  density  of 
shot  noise  was  reported  by  Richter  and  Smits  [20]  in  1974. 
Their  approach  was  to  approximate  the  characteristic  func¬ 
tion  by  piecewise  polynomial  segments,  which  could  then 
be  inverse  Fourier  transformed  in  closed  form.  In  1975  Fos- 
chini  ei  al.  [9]  gave  a  detailed  analysis  of  the  detection  of  a 
shot-noise  signal  in  the  presence  of  additive  white  Gaussian 
noise.  They  employed  several  approximations  in  order  to 
obtain  manageable  expressions  for  the  likelihood  function, 
from  which  they  could  gain  insight  into  the  general  struc¬ 
ture  of  the  optimum  detector.  The  1976  paper  of  Mazo 
and  Salz  [16]  analyzed  the  performance  of  integrate-and- 
dump  filters.  They  also  obtained  exact  formulas  for  the 
probability  distribution  of  the  {/!„}  in  physical  avalanche 
diodes.  Further  work  on  computing  the  density  of  shot 
noise  appeared  in  the  1978  paper  by  Yue  ei  al.  [23].  Their 
approach  was  to  approximate  the  tails  of  shot-noise  densi¬ 
ties  by  a  weighted  sum  of  normal  densities.  In  1984  Morris 
[17]  considered  an  imaging  model  in  which  the  photon  lo¬ 
cations  {7j,}  were  observed  and  passed  through  a  linear 


filter  matched  to  the  underlying  Poisson-process  intensity 
A( •).  He  then  studied  hypothesis  testing  based  on  the  re¬ 
sulting  shot-noise  random  variable.  In  1988  Kadota  [14] 
reported  approximately  optimum  detection  of  determinis¬ 
tic  signals  in  Gaussian  and  compound  Poisson  noise.  His 
model  included  a  noise  process  that  consisted  of  samples 
of  a  Gaussian  process  where  the  sample  times  were  Pois¬ 
son  distributed.  In  1990.  Lowen  and  Teich  [15]  considered 
shot-noise  processes  with  a  power-law  impulse  response  to 
obtain  1//  noise.  They  assumed  that  the  underlying  Pois¬ 
son  process  had  a  constant  intensity,  and  they  showed  that 
such  shot-noise  distributions  need  not  converge  to  a  normal 
law  as  the  intensity  increases.  In  1991  Hero  [13]  approxi¬ 
mated  the  likelihood  ratio  for  observing  a  shot-noise  pro¬ 
cess  in  the  presence  of  additive  white  Gaussian  noise.  His 
approximation  became  more  accurate  as  impulse  response 
of  the  shot  noise  became  more  narrow  and  more  closely 
resembled  the  underlying  point  process.  Most  recently, 
Heist  rom  and  Ho  [12]  have  successfully  applied  numeri¬ 
cal  contour  integration  to  the  inversion  of  the  shot-noise 
characteristic  function  to  obtain  the  shot-noise  cumulative 
distribution  in  the  case  when  additive  Gaussian  noise  is 
present. 


IV  Shot-Noise  Characteristic  Functions 

In  the  remainder  of  the  paper  we  fix  t.  set  g(r)  :  = 
h(t.  r).  and  focus  on  the  random  variable 

V  :=  A„g(Tl/). 

U 

The  characteristic  function  of  V  is  [22,  p.  170], 


As  suggested  by  the  dr  in  (10),  we  are  implicitly  as¬ 
suming  that  the  underlying  Poisson  process  lives  on  some 
finite-dimensional  euclidean  space.  For  temporal  Poisson 
processes,  the  integral  in  (10)  is  over  IR  =  (-00,00).  For 
twodimensional  spatial  Poisson  processes,  we  have  a  dou¬ 
ble  integral  over  the  euclidean  plane  IR2.  In  any  case,  using 
Fubini’s  Theorem,  we  can  rewrite  (10)  as 


t/'M  =  A(r)[e>“'»(™*'  -  1  }dr 

assuming  that 


i  J  A(r)|ej“j(™‘'  -  1 1  dr 


<  oc. 


(11) 


(12) 


In  order  that  (12)  be  true  even  when  ,4„  is  a  constant 
random  variable,  we  assume  that 


/A(r)|e— -l|dr  <  00.  for  all  u  €  IR.  (13) 

Remark:  A  sufficient  condition  for  (13)  to  hold  is  that 
the  integral  of  A  over  the  support  of  g  be  finite.  For  exam¬ 
ple,  consider  a  temporal  Poisson  process  with  A(r)  =  0  for 
t  <  0.  Let  h  be  a  causal  impulse  response,  i.e.,  h(t,  r)  =  0 
for  r  >  t.  Then  with  g(r)  =  h(t,  r)  for  some  fixed  t  >  0, 
the  integral  in  (13)  becomes 

A(r)|eru,/,((,r)  _  j|dr 

Clearly,  this  integral  will  be  finite  for  every  t  >  0  and  every 
^  €  IR  assuming  only  that  A  is  locally  integrable. 

In  view  of  (1 1),  we  set 


fM  :=  E[e-'“v]  =  exp(t£>(w)), 

where 

vM  :=  Jx(r)M«g(T))-l]dr,  (10) 
and  r'a  is  the  common  characteristic  function  of  the  {A„}. 

Remark:  If  we  had  added  Gaussian  noise  in  the  defi¬ 
nition  of  V,  then  the  characteristic  function  of  Y  would 
have  decayed  like  e~°  w  ^ ,  where  <r2  is  the  variance  of  the 
Gaussian  noise.  In  this  case.  Y  would  have  an  infinitely 
differentiable  density,  and  the  method  of  Section  II  would 
apply 

The  purpose  of  this  section  is  forgive  conditions  under 
which  we  can  write  v(u/)  =  —  B  +  tr(wi),  where  — *  0 

as  M  —  x.  It  will  then  follow  that  the  density  of  V 
contains  an  impulse  of  strength  e~B  at  the  origin.  We 
will  also  identify  the  constant  B  and  the  inverse  Fourier 
transform  of  v. 


t/’0(w)  :=  J  A(r)[e’u,,'T*  —  l]dr  (14) 

so  that  we  can  write 

i>(u)  -  E[t/i0(wA1/)].  (15) 

To  analyze  (14),  we  proceed  as  follows.  YVe  first  define  two 
measures, 

A (D)  :=  (  A (r)dr, 

Jd 

where  D  is  any  Borel  subset  of  the  euclidean  space  on 
which  the  underlying  Poisson  process  lives,  and 

F(C)  :=  A({r  :  g(-r)  £  C}), 

where  C  is  any  Borel  subset  of  IR. 

Remark:  If  A(  D)  <  oc  for  some  set  D,  then 

A({r:g(r)€C}n£>) 

A  (D) 


m 


can  be  regarded  as  the  conditional  probability  of  the  event 
{^(jHi)  €  C]  given  that  exactly  one  Poisson-distributed 
point  occurs  in  D.  If  {r:  g(r)  £  C}  C  D ,  then  ( 16)  reduces 
to  T(C)/A(D). 

We  now  apply  the  change-of-variable  formula  for  mea¬ 
sures  [4,  p.  219,  Theorem  16.12]  to  (14)  to  obtain 

if-oM  =  J[eju9{T)-l}dA(r) 

=  f°°  [ei“B  -  l]rfr(0). 

J  —  oo 

Since  the  integrand  is  zero  for  9  =  0,  we  can  write 

/CM 

[e1*9  —  1]  dro(0), 

-OO 

where  To  is  the  measure  defined  by 

r0(C)  :=  r(c\{o}). 

If  To  is  a  finite  measure,  i.e.,  if 

r0(IR)  =  /  A(r)  dr  <  oo, 

J{r  j(r)*0} 

then  we  can  set  Bo  '■=  To(IR)  and  write 

/oo 

e]wtd  r0(fl). 

-OO 


If  r0  has  a  density  70.  then  we  can  write 

V-o(w)  =  -B0+  j  e^'%  (0)d6. 

J  —  CO 


(17) 


Proposition  2:  If  To  is  a  finite  measure  with  density  70, 
then 

t(^)  =  -B+  l  e9*9-,  (0)dO, 

J  —  OO 


where 

and 


B  :=  B0  P(AV  >0), 
7o  (9/ Ay). 


7 (*)  :=  E 


LM»>o} 


(18) 


If  we  now  set 

£(->)  :=  B  +  xi'(uj)  =  f  eJU,9y(9)d9, 

J  —  OO 


then  ip(uj)  — *  0  as  |w|  — *  00  by  the  Riemann-Lebesgue 
Lemma.  It  then  follows  that  <p(lj)  — -  e~B  as  |w|  —  00,  and 
therefore  the  density  of  Y  contains  an  impulse  of  strength 
e~B  at  the  origin. 

An  analogous  development  can  be  carried  out  for  the 
moment-generating  function  of  Y, 

m(s)  :=  E[elK]  =  exp(F(s)), 

where  for  s  6  C,  k(s)  :  =  E[£a(sA„)], 

Ms)  ==  J  A(r)[e’*<T'-l]dr, 

and  we  assume  that  for  all  a  €  Bt, 


E  [  A(t -)e°^A-  dr 

,J{t  g(T)?0} 


<  OO. 


Note  that  when  a  =  0,  this  reduces  to  the  weaker  assump¬ 
tion  that  r0(!R)  <  00. 

A.  Examples 

Let  ||  - 1|  be  a  norm  on  IR2.  Fix  any  r  >  0,  and  let  q  map 
[0,  r]  onto  [0, 1],  where  q  is  nonnegative,  continuous,  and 
strictly  decreasing  so  that  ':[(),  1]  —  [0,r]  exists  and 
is  nonnegative,  continuous,  and  strictly  decreasing.  Set 
g(r)  =  ?(IMI)  for  IMI  <  r  and  set  g(r)  =  0  for  ||r||  >  r. 
Since 

g(T)  >  9  <=>  ||r||  <  q~l(9),  0  <  0  <  1, 

the  definition  of  To  yields 

ro(0,oo)  =  f  A (r)dr,  0  <  9  <  1, 

JD(B) 


where 


D{6)  :=  {r€W  :\\r\\<q-l(e)}. 


is  integrable;  the  symbol  1  is  the  indicator  function  of  the 
specified  set. 

Proof:  First,  since  70  is  integrable,  Tonelli’s  Theorem 
and  a  change  of  variable  imply  that  7  is  integrable.  Next, 
substitute  ( 17)  into  ( 15)  and  break  up  the  expectation  over 
the  events  (.4,  =  0}  and  {.4„  >  0}.  Make  a  change  of 
variable  and  then  use  Fubini’s  Theorem.  □ 


(19) 


Since  g  is  nonnegative,  Fo(0,oo)  =  Bo  for  9  <  0,  and  since 
the  maximum  value  of  g  is  1,  Fo(0, 00)  =  0  for  9  >  1.  The 
density  of  To  is  7o(0)  =  —  (d/d0)r<)(0, 00).  We  next  let 

Hr)  =  u(||r||2), 

where  u:[0,oo)  — *  [0, 00).  If  ||  ■  ||  is  the  usual  euclidean 
norm,  then  changing  to  polar  coordinates  yields 

/■«-'(«> 

To(0,oo)  =  2x  /  u(p2)pdp,  0  <  9  <  1. 

Jo 

If  U  is  an  antiderivative  of  u,  then 


—  f/(p2)  =  2pu(p2), 


7cZ 


and 


ro(0,oo)  =  ir[^(q-‘(^)2) -t/(0)],  O<0<1. 

For  example,  the  Gaussian  intensity,  A(r)  =  exp(-||r||2), 
corresponds  to  u{p)  =  e~p  and  U(p)  =  —t~p.  In  this  case, 


To(0,  oc)  =  tt(  1  —  e 


0  <  9  <  1. 


q(x)  = 


(20) 


rl{9) 

and  if  A(  r)  =  1 

rn(0.  *)  = 

and 

7o(fl)  =  \ 


j  1  -  sJFfb,  0  <  e  <  q(z), 

\  v/(l  -9)/a,  q(:)<9<  1, 

*{l-2y/0/b  +  0/b),  0  <9<q(z), 

7r(  1  —  6)/a,  q(z)  <  9  <  1, 


t((^/6)-'/2-  l)/6,  0  <e<q(z), 

tt /a,  q(z)  <  6  <  1. 

I'sing  ( 17),  we  find  that  for  ui  >  0, 

e;w  _ 


Vo(w)  —  w 


+ 


juia 

1  _  * 1 


j*jb 


+ 


Wt1-1}’  <*> 


where  r  =  \f‘2q(z)<^ j t,  and 

Fr(x)  :=  f  exp (j$02)d0. 
Jo 


The  real  and  imaginary  parts  of  Fr(x)  are  Fresnel  integrals, 
and  are  easily  computed  with  the  NAG  library  subroutines 
S20ADF  and  S20ACF,  respectively.  In  addition,  for  s  >  0, 
we  can  write 


*o(s) 


= 


The  Cauchy  intensity,  A(r)  =  l/(  1  +  ||r||2),  corresponds  to 
u(p)  =  1/(1  +  p)  and  U (p)  =  ln(  1  +  p),  yielding 

ro(0,oc)  =  jt  ln(  1  +  q~'(0)2),  O<0<1. 

The  constant  intensity  A(r)  =  1,  corresponds  to  u(p)  =  1 
and  U(p)  =  p ,  resulting  in 

r0(fl.oo)  =  tt?-1^)2,  0  <  6  <  1. 

Example  2:  Consider  the  piecewise  quadratic  function 

1  —  ax2,  0  <  x  <  z, 

6(1  —  x)2,  z  <  x  <  1, 


e*  - 
sa 
1  _ 

+  7b 


where 


Daw(x) 


+  2_y^"Daw( \/sq(z) )  -  1  j. 

=■  •-/* 


e9*  d9 


is  Dawson’s  integral,  and  is  easily  computed  with  the  NAG 
library  subroutine  S15AFF. 

Example  3:  Let  q  be  as  in  (20)  with  c  =  1  and  a  =  1. 
This  results  in 


where  a  >  0,  6  >  0,  and  0  <  z  <  1  are  given.  If  0  <  r  <  1, 
we  see  from  (20)  that  q( 0)  =  1,  <j(l)  =  0,  and  <j'(0)  = 
q'(  1)  =  0.  If  a  =  6/(6  —  1)  and  r  =  b/(a  +  6),  then  q  will 
be  continuously  differentiable  at  z,  and  thus  q  will  be  a 
quadratic  spline  on  [0,  1  J,  an  example  of  which  is  shown  in 
Fig  4  Now, 


q-'( 9)  =  v/1  -  0,  0  <  9  <  1. 

If  A(r)  is  the  Cauchy  intensity  mentioned  earlier,  then 
ro(0,oc)  =  7rln(2-0),  0  <  6  <  1, 

70 (9)  =  »/(2-0),  0  <  6>  <  1, 

and  for  uj  >  0, 


tl’o(uj)  =  — 7r  In  2  + 


f1  e]w 

*  Jo  2- 


=  —  rr  In  2  Tire7 


2ui 


d9 

7w 


L  — d 


The  real  and  imaginary  parts  of  this  last  integral  are  given 
by 

Ci(2ui)  -  Ci(w)  and  Si(w)  -  Si(2w), 
respectively,  where 

Ci(w)  :=  C  +  lnw  +  /  COS0-  1  d9, 

Jo  9 

C  a  0.5772  is  Euler’s  constant,  and 

.  [w  sin  9 

Si(w)  :=  J f  —Q~d6' 

The  functions  Ci  and  Si  are  easily  computed  with  the  NAG 
library  subroutines  S13ACF  and  S13ADF,  respectively.  In 
addition,  for  s  >  0, 


k0(s)  =  -t In 2  +  tre2'[E](s)  -  £'i(2s)], 


where 


Ex(x) 


f°°  e~ 9 

-  I  T* 


x  >  0, 


is  the  exponential  integral,  and  is  easily  computed  with  the 
NAG  library  subroutine  S13AAF. 


is  nonatomic,  and 


g(r)  =  g(x 


■-  ■  { 


Example  4  (Moms  [17]/  For  r  =  ( x,y )  G  IR2,  set 

-  1*1.  1*1  <  1.  kl  <  1, 

otherwise, 

and  let  A(r)  =  (Nr/2)g(r).  Then 

r0(fl,oo)  =  AMI -0)(1 +0),  0  <  <  1, 

and 

ToW  =  2  Nr0,  0  <  0  <  1. 

Note  that  B0  =  r0(IR)  =  Fo(0,  oo)  =  Nr-  Using  (17), 

rJU,/2  /  ,sin(w/2) 


t/'oM 


Bn 


-1  + 


w/2 

-;'cos(w/2) 

Also,  for  s  /  0, 

*o(«)  =  So 


sin(w/2)  +  j 


w/2 


w  /  0. 


e*  eJ  -  1 

-  I  +  2[  — =— 


q-'(0)  = 


-ln0 


,  0  <  0  <  1. 


If  A(r)  is  the  Gaussian  intensity  mentioned  earlier,  then 
ro(0.oo)  =  tt(1  -01/p),  0  <  0  <  1, 


<$(w)  =  f  eJU'z  dfi(x) 

J  —  OO 

=  E  [e^l^o,] 

=  v?M-e-B 
=  e-*[e*(«)  _  !]. 

Since  ^>(w)  — *  0,  $(w)  s=  e-Bt/>(w)  for  large  |w|.  Further¬ 
more,  \p  can  decay  very  slowly  for  shot-noise  processes. 
Recalling  Rummer’s  Comparison  Method  [11,  p.  195],  we 
write 

*(w)  =  e-B[e^> -^(w)-l]  +  e-B^(w), 

where  the  term  in  brackets  decays  like  V'(w)2/ 2.  Substi¬ 
tuting  ‘his  into  (3),  we  find  that 


H(y,oc)  =  e  B[Gc{y)  +  r/(y)], 


(23) 


Example  5:  Let  q{x)  =  exp (-pr2),  where  p  >  0.  If  we 
regard  q  as  mapping  [0,oo)  onto  (0,  1],  then 


where 


Gc(y)  :=  lim  V  bn[e^nr/L) -^(nir/ L)- \]e~>n"'lL, 

L  —  oo  ^  * 
n  =  —  oo 

and 

OO 

rj(y)  :=  lim  V]  bnrp(me/ L)e~jni,ylL . 

L—*oo 


and 

7o(»)  =  ~0’~\  0  <  0  <  1. 

P 

Clearly,  if  p  =  1/n  for  some  positive  integer  n,  then  on 
[0.  1],  7U(0)  is  a  polynomial  of  degree  n  —  1,  and  (17)  can  be 
evaluated  in  closed  form  using  integration  by  parts.  How¬ 
ever,  if  p  >  1,  then  7o(0+)  will  not  be  finite;  for  example, 
if  p  =  2,  7o  will  exhibit  1  j \f0  behavior  as  in  Example  2. 

Remark  If  z  =  1  and  a  =  1  in  Example  2,  and  if 
p  =  1  in  Example  5,  then  the  resulting  70  functions  are 
the  same,  7o(0)  =  t.  If  p  —  1/2  in  Example  5,  and  if 
,Vr  =  it  in  Example  4,  then  the  70  functions  are  again  the 
same,  7o(0)  =  2*0. 


V.  Shot-Noise  Cumulative  Distributions  and 
Densities 

Let  F  denote  the  cumulative  probability  distribution  of 
Y.  We  assume  that  F  is  continuous  everywhere  except  the 
origin,  where  it  has  a  jump  discontinuity  of  size  e_B.  Then 
the  measure 


p(C)  :=  P(V  G  C  and  Y  f  0)  (22) 


From  (22)  and  (23)  it  follows  that 


P (Y  >  y) 


«  B[Gc(y)  +  R(y)],  y>  0, 

e-B[Gc(y)  +  rj(y)  +  1],  y  <  0. 


Fortunately,  we  can  evaluate  the  limit  defining  rj  to  obtain 


7  (0)d0-, 


this  follows  by  applying  (3)  to  the  measure  whose  density 
is  7(0)  and  by  noting  that  the  corresponding  characteristic 
function  is  given  by  (19). 

We  now  observe  that  if  Gc  is  absolutely  continuous,  then 
F(y)  =  1  —  P(y  >  y)  has  density 


/(y)  =  e 


- 


Hy)  +  7 (y)  -  j-Gc(y) 

dy 


(24) 


Remark:  Equation  (24)  generalizes  to 

m—  1 


/(y)  =  C 


- 


Hy)  +  7(y)  +  ^2  7 *’(y)  -  ^(y) 


i  =  l 


where  7“  is  the  convolution  of  7  with  itself  i  times,  and 

00  m  ~ 

Gcm(y)  ■  =  lim  T  6n[e*‘""t>-£  +l2E(LL]e-i"*/L. 

*=° 

To  approximate  the  last  term  in  (24),  we  fit  a  cubic 
spline  to 

.v 

Gci  .v(y)  :=  Y1  K[e*(n'IL)  -iinxlL)-\}e->n">lL 

n = -  S 

(25) 

and  then  differentiate.  Note  that  even  if  V’(^)  decays  like 
l/v-CJ,  as  it  does  in  Example  2  if  we  assume  Av  —  1  with 
probability  1,  the  coefficients  in  (25)  decay  like  1/n2,  and 
so  as  iV  —  00,  Gcl  n  converges  to  a  continuous  function. 
Thus,  for  finite  A\  (25)  will  not  exhibit  Gibbs  phenomenon. 
Of  course,  if  we  were  to  differentiate  (25),  the  new  series 
coefficients  would  decay  like  1/n,  and  Gibbs  phenomenon 
would  appear.  By  fitting  a  spline  with  appropriately  cho¬ 
sen  knots  to  Gci  v  and  then  differentiating,  we  can  avoid 
the  Gibbs  oscillations. 


A.  Error  Bounds 

By  taking  L  finite  in  the  first  sum  in  (23)  we  see 
that  there  are  two  sources  of  error.  The  first  one  can 
be  bounded  by  looking  at  the  measure  corresponding  to 
^(w)  -  eB,  i.e.,  to  /i  itself  as  was  done  in  Section  II.  The 
second  source  can  be  bounded  by  looking  at  the  measure 
corresponding  to  ip(u),  i.e.,  to  the  measure  whose  den¬ 
sity  is  7.  Now,  if  |y(r)|  is  bounded  by  some  constant, 
say  ymax,  then  7 o(0)  =  0  for  |0|  >  ymax.  Hence,  if  is 
a  bounded  random  variable,  say  Au  <  a,  then  by  (18), 
7(0)  =  0  for  \6\  >  aymax.  So,  if  L  +  L?  >  a ymax,  and  if 
—  L  +  L\  <  — oymax ,  there  will  be  no  error  due  to  the  second 
source.  This  will  be  the  case  in  the  following  examples. 


B  Xumertcal  Examples 


In  the  following  examples  we  take  Au  =  1  so  that  = 
Co,  k  =  k0,  and  7  =  70. 

Example  2  Continued:  Let  a  =  b  =  2  and  ;  =  1/2.  Then 
y(r)  =  y(||r||),  where  q  is  the  quadratic  spline  shown  in 
Fig.  4.  We  first  consider  the  error  bounds.  Since  g(r)  >  0, 
/i(—  oc,  0)  =  0,  and  so  the  upper  bound  in  (4)  is  zero.  We 
now  take  L  =  7,  L\  =0,  and  Ln  =  5.  To  analyze  the  lower 
bound  in  (5),  write 


H(L  +  L  1.3c)  = 


< 


fi(L.oc) 

P( Y  >  L,Y  £  0) 

P(V  >  L),  since  L  >  0, 

e~,Lm(s) 

e-[sL-k0it)] 


It  is  easily  verified  numerically  that  the  maximum  value  of 
7s  -  k0(s)  >s  greater  than  12  and  occurs  at  s  ss  3.  This 
yields  the  bound  e-12  <  6.2  x  10-6.  Turning  now  to  (25), 
we  take  N  =  256.  An  FFT  provided  93  samples  of  GCL  N  in 
the  interval  [0,5.01].  We  then  fit  aspline  with  64  uniformly- 
spaced  knots  in  the  proper  subinterval  [.01,5].  A  plot  of 
minus  the  derivative  of  this  spline  is  shown  in  Fig.  5.  There 
appear  to  be  corners  at  1/2  and  1.  To  confirm  this,  we 
plotted  minus  the  second  derivative  in  Fig.  6.  Clearly, 
there  is  a  large  downward  jump  at  1  but  only  a  corner  at 
1/2.  In  addition,  although  it  is  difficult  to  see,  there  is  a 
jump  at  2;  just  as  there  are  oscillations  on  both  sides  of 
the  jump  at  1,  there  are  telltale  oscillations  on  both  sides 
of  2.  If  we  add  a  single  knot  at  .5  and  double  knots  at 
1  and  2,  and  recompute  the  curve  in  Fig.  6,  we  obtain 
the  curve  in  Fig.  7.  Notice  how  the  oscillations  around 
1  and  2  have  disappeared.  Based  on  these  observations, 
we  selected  a  new  set  of  21  nonuniformly  spaced  knots, 
including  a  single  knot  at  .5  and  double  knots  at  1  and 
2.  Using  the  original  93  samples  of  GCL  v,  we  fitted  a  new 
spline  with  the  set  of  21  knots.  After  differentiating  and 
negating,  we  obtained  the  curve  shown  in  Fig.  8.  If  we  add 
7(y)  to  this  curve  and  multiply  by  e~B  (cf.  (24)),  we  obtain 
the  graph  in  Fig.  9.  For  comparison,  we  have  also  plotted 
with  a  dashed  line  the  256-point  FFT  approximation  of 
f{y)  obtained  via  equation  (8).  The  impulse  at  the  origin 
is  not  shown  in  the  figures. 

Example  3  Continued:  Proceeding  as  in  the  previous  ex¬ 
ample,  again  with  L  =  7,  L\  =  0,  and  L2  =  5,  we  ob¬ 
tain  the  error  bound  e~8  65  <  1.8  x  10~4.  With  N  =  256 
there  were  93  samples  of  G\  N  in  the  interval  [0,5.01].  We 
started  by  fitting  a  spline  with  64  uniformly  spaced  knots 
in  the  subinterval  [  01,5].  Upon  inspection  of  the  result 
(not  shown),  we  added  double  knots  at  1  and  2  for  a  total 
of  68  knots.  A  plot  of  minus  the  derivative  of  this  spline 
is  shown  in  Fig.  10.  Notice  that  this  curve  is  continuous 
with  corners  at  1  and  2.  Based  on  this  curve,  we  selected 
17  nonuniformly  spaced  knots,  including  double  knots  at 
1  and  2.  Using  the  original  93  samples  and  the  new  set 
of  17  knots,  we  generated  a  new  spline  approximation  of 
G\  N.  A  plot  of  e~B  times  the  quantity  7(y)  minus  the 
derivative  of  the  17-knot  spline  is  shown  in  Fig.  11.  For 
comparison,  we  have  again  plotted  with  a  dashed  line  the 
corresponding  256-point  FFT  approximation  of  f(y). 

Example  4  Continued:  Let  Nr  —  2.  With  L  =  8,  Li  =  0, 
and  L2  =  6,  we  obtained  an  error  bound  of  e-9  5  <  7.5  x 
10-5.  With  N  =  256,  there  were  98  samples  of  G\  N  in 
the  interval  [0,6.06].  We  started  by  fitting  a  spline  with 
50  uniformly  spaced  knots  in  the  subinterval  [  01,6].  Upon 
observing  the  result  (not  shown),  we  generated  a  set  of  16 
nonuniformly  spaced  knots,  including  a  double  knot  at  2. 
A  plot  of  e~B  times  the  quantity  7 (y)  minus  the  derivative 
of  the  16-knot  spline  is  shown  in  Fig.  12.  The  jump  at  1 


is  due  to  the  fact  that  7  is  not  continuous  there.  Again, 
for  comparison,  we  have  plotted  with  a  dashed  line  the 
corresponding  256-point  FFT  approximation  of  f(y).  In 
[17,  Fig.  4(a)],  Morris  obtained  a  similar  solid  curve  using 
an  8192-point  FFT. 


VI.  Approximation  of  70  and  Vo 

As  we  have  seen,  in  order  to  employ  the  methods  of 
the  previous  section,  we  must  be  able  to  compute  both 
7  and  4>,  which  can  be  expressed  in  terms  of  70  and  V 0, 
respectively  (cf.  (18),  (19)  and  (15)).  In  this  section  we 
discuss  an  approximation  scheme  to  compute  both  70  and 
Vo- 

We  first  consider  the  Fourier  integral  in  (17).  Without 
loss  of  generality,  we  restrict  attention  to  the  integral  over 
the  right  half  line, 

e^7o  (0)d9.  (26) 

VV'hile  70  is  integrable  since  it  is  the  density  of  a  finite 
measure,  it  is  still  possible  to  have  7o(0+)  =  00  as  we 
saw  in  Example  2.  Thus,  it  may  be  difficult  to  compute 
this  integral  numerically.  We  therefore  propose  the  appli¬ 
cation  of  integration  by  parts.  However,  before  doing  so, 
it  is  convenient  to  introduce  the  following  notation.  For 
9  >  0,  let  Ho{9)  .=  Told,  00).  (For  9  <  0  we  would  use 
H[)(9)  =  ro(-oo,0).)  Since  To  is  a  finite  measure  with 
density  70,  it  immediately  follows  that  for  6  >  0,  Ho  is 
bounded,  continuous,  nonincreasing,  and  decays  to  zero  as 
9  goes  to  infinity.  Also,  Hq(9)  =  -70 (9)  for  almost  every 
9  >  0.  Now  set 

Hn+l(0)  :=  -J  //„(<) d<,  9  >  0.  (27) 

Since  Ho  is  nonnegative,  Hn  is  always  of  one  sign. 

Proposition  3:  Let  G+  :  =  {r  :  g(r)  >  0}.  For  v  >  0,  set 
Sa(9.  v )  :  =  l[o, »)(#),  and  for  n  >  1,  set 

(9  —  v)n/n!,  O<0<t>, 

0,  9  >  v. 


Then 

fln(9 )  =  /  \(r)S„(9,  g(r))  dr,  9>  0.  (28) 

Jg+ 

Furthermore,  Hn(9)  is  finite  for  9  >  0  if  and  o;  ly  if  Hn( 0) 
is  finite;  i.e.,  if  and  only  if 


< 


00. 


If  l!n( 0)  is  finite,  then  Hn  is  continuous,  Hn(9)  —  0  as 
0  — *  o c,  and  H'n(9 )  =  Hn-\(9)  is  integrable. 


Proof:  Use  induction  and  Toneili’s  Theorem.  □ 

Remarks:  (»)  Since  we  are  assuming  that  r0  is  a  finite 
measure,  /G+  A (r)dr  <  00.  Thus,  if  g  is  bounded,  Hn( 0) 
is  finite  for  all  n. 

(it)  For  fixed  9,  So(9,v)  is  a  discontinuous  function  of 
t>,  Si(0,v)  is  a  continuous  function  of  v,  and  for  n  >  2, 
S„(0,  t>)  is  an  n  —  1  times  continuously  differentiable  func¬ 
tion  of  v. 

Corollary  4:  If  Hn+ 1(0)  is  finite,  then  the  integral  in  (26) 
is  equal  to 

n  roo 

£(-jw)*//fc(0)-(-jW)"+1  /  ei“eHn(6)d9,  (29) 

i  =  0  *'° 

where  Hn  is  n  times  continuously  differentiable. 

Rather  than  compute  the  integral  (26)  numerically,  it 
may  be  advantageous  to  use  (29)  instead,  since  Hn  is 
smoother  than  70.  We  propose  using  H\  in  the  follow¬ 
ing  approximation  scheme  for  computing  both  70  and  Vo- 
Assume  //2(0)  is  finite.  We  compute  H\(0)  at  just  enough 
points  9  so  that  we  can  obtain  a  good  cubic  spline  approx¬ 
imation  of  H\.  Denote  this  approximation  by  H\.  Then 
—H"  will  provide  a  piecewise  linear  approximation  of  70. 
Of  course,  the  approximation  will  break  down  near  the  ori¬ 
gin  if  7o(04-)  =  00.  To  approximate  ^(w)  atu>  =  nir/L  for 
n  =  0, . . . ,  JV  —  1,  we  use  (29)  as  follows.  Set  AS  =  2 L/M 
for  some  M  >  IV.  Let  Hi  denote  the  piecewise  linear  ap¬ 
proximation  of  Hi  satisfying  Hi(kA9)  =  Hi(kA0).  Then 
if  we  substitute  Hi  into  (29)  and  use  integration  by  parts, 
we  obtain 

M-l 

H o(0)  +  ct[rJU'(t+1>A<’  - 

k= 0 

where  ct  is  the  slope  of  Hi  on  the  interval  ( k  A 9,  (k  + 
1)A0).  Assuming  M  is  large  enough  that  Hi((M—l)A9)  = 
Hi(M  A 9)  —  0,  we  will  have  cat-  1  =  0.  Then  since  AS  = 
2 L/M,  if  u)  =  nir/L,  the  summation  can  be  evaluated  for 
n  =  0, . . . ,  M  -  1  with  a  pair  of  A/-point  FFTs. 

In  the  continuation  of  Example  2,  we  computed  Vo  ex¬ 
actly  using  (21).  We  have  also  repeated  those  calculations 
replacing  the  true  values  of  \po(nrr/ L)  by  the  just-described 
approximation.  The  results  were  graphically  indistinguish¬ 
able  from  the  solid  curves  already  shown.  The  particulars 
were  as  follows.  Using  64  uniformly  spaced  knots  on  the 
subinterval  [10— 7 , 1  —  10-7],  we  fitted  a  cubic  spline  to  96 
uniformly  spaced  samples  of  Hi  from  the  interval  [0, 1].  To 
approximate  used  M  =  1024.  Similar  results  were 

also  obtained  using  only  46  nonuniformly  spaced  samples 
of  Hi  and  16  nonuniformly  spaced  knots.  (Of  course,  on 
a  fine  enough  scale,  the  approximations  of  70  break  down 
near  the  origin.  For  these  particular  approximations,  the 
breakdown  is  not  visible  graphically  for  9  >  .05.) 


Appendix 

Proof  of  Lemma  1 


References 


Our  approach  is  a  generalization  of  [1].  First  note  that 
A<(0, oo )  =  l(0,oo)(*)<M*)-  Now,  for  any  L  >  0,  de¬ 

fine  the  periodic  pulse  train  3l  with  period  2 L  by  specify¬ 
ing  its  values  on  (  —  L,  L\  to  be 


3l(x) 


1,  0  <  x  <  L, 

0,  -L<x<  0. 


Since  — ■  l(o,cc)(J)  for  all  x  as  L  —  oo,  the  domi¬ 

nated  convergence  theorem  implies 


3l(x)  d/i(x)  —  /i( 0,  oc). 


(30) 


Below  we  show  that  3l(x)  dfi(x)  is  given  by  the  in¬ 
finite  series  in  (1).  Before  doing  so,  we  first  give  another 
proof  of  (30)  that  yields  a  bound  on  the  error  when  L  is 
finite.  Observe  that  for  all  x. 


3l(*)  <  1(0, cci(-r)  +  !(-«, _/.](*) 


and 


■h(-r)  >  1(0. L^)  =  l(0.rc)('r)  —  l(L.oo)(r)- 

Hence,  the  difference, 


3L(x)d/j(x)  -  /t(0.  x), 


is  upper  bounded  by  /i(  — x,  -/,]  and  lower  bounded  by 
-//(L.x).  Since  ft  is  a  finite  measure,  ^(L.oc)  and 
/<(-x.  - L ]  both  converge  to  zero  as  L  goes  to  infinity. 

The  next  step  in  the  proof  is  to  approximate  3l  by  its 
Fourier  series.  Set 


v 

3l,n(*)  :=  Y,  6nejn”/i, 

n  —  -S 


where  the  bn  are  given  by  (2).  Now,  as  N  —  x,  3iy(x)  — 
Ji  ( x )  for  all  x  except  at  points  of  discontinuity  of  3l  Since 
/i  is  finite  and  nonatomic,  and  since  3i,y(x)  is  uniformly 
bounded  in  ,V  and  x  (with  L  fixed)  [8,  p.  151,  eq.  (10.1.5)], 
the  dominated  convergence  theorem  yields 


3f.(x)  dfi(x) 


e]nrz/L  dn(x) 


rc 

Y  bnQ>{nir/L).  □ 

n~  —  ~c- 
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I.  Introduction 


When  images  are  produced  under  high-level  illumination,  they  are  often  modeled  as  the  sum  of 
a  “signal”  image  plus  signal-independent  Gaussian  noise.  However,  when  an  image  is  formed  under 
low-light-level  illumination,  it  can  be  better  modeled  as  a  filtered  point  process,  also  known  as  shot 
noise,  described  as  follows. 

The  process  of  image  acquisition  consists  of  the  measurement  of  arriving  photons  in  the  image 
plane.  In  practical  situations  it  is  difficult  to  measure  the  exact  location  of  the  photons  since 
these  are  filtered  by  the  finite  response  of  the  imaging  device.  What  is  measured,  instead,  is  the 
superposition  of  the  responses  of  the  imaging  system  to  each  arriving  photon.  This  superposition, 
measured  at  a  point  x  6  1R2,  can  be  described  by  the  random  variable 

Z(x)  =  £h(x-Xu),  (1.1) 

where  { X „}  denotes  the  set  of  positions  at  which  photons  are  detected,  and  h  represents  the  impulse 
response  or  point  spread  function  of  the  imaging  device.  A  block  diagram  of  the  imaging  system  is 
shown  in  Fig.  1. 

We  consider  the  following  hypothesis-testing  problem.  Under  hypothesis  H,,  i  =  0,1,  the  {AV} 
are  points  of  a  two-dimensional  Poisson  process  with  nonnegative  intensity  A;(x),  x  €  1R2. 

Special  cases  of  this  problem  that  have  previously  been  addressed  include: 

1.  The  exact  locations  of  the  photons  are  available.  This  is  equivalent  to  the  point  spread 
function  h(x)  being  equal  to  a  Dirac  delta  function  6(x),  in  which  case  one  has  a  detection 
problem  with  Poisson-process  observations.  When  one  has  Poisson-process  observations,  the 
likelihood  ratio  (LR)  is  of  course  well  known  [9,  p.  94],  In  [4],  a  suboptimal  detection  scheme 
that  could  be  easily  implemented  was  considered.  This  led  to  the  consideration  of  a  correlation 
detector  in  which  the  Poisson- process  observations  were  passed  through  a  linear  filter  (taken 
to  be  one  of  the  intensity  functions)  and  sampled.  This  led  to  a  hypothesis  test  based  on  a 
single  shot-noise  random  variable. 

2.  Counts  of  photons  in  disjoint  regions  are  available,  and  hence  one  is  faced  with  a  Poisson 
counting  process  detection  problem.  This  case  can  be  regarded  as  a  filtered  Poisson  process 
with  special  form  of  h,  followed  by  sampling.  The  LR  is  well  known  for  this  special  case  also  [9, 
p.  94].  In  [10],  a  correlation  scheme  was  used  for  classification.  The  Poisson  counting  process 
observations  (with  counts  being  either  0  or  1  with  high  probability)  were  cross-correlated 
with  various  reference  functions.  Three  reference  functions  were  considered,  one  of  which  was 
constructed  so  that  the  value  of  the  correlation  between  that  function  and  the  observed  image 
approximates  the  value  of  the  logarithm  of  the  LR. 

Little  work  has  been  done  for  the  more  general  case  of  filtered  point-process  observations  due  to 
the  difficulty  of  computing  the  density  and  distribution  functions  involved  [3]. 

II.  Hypothesis  Testing  with  Shot-Noise  Observations 
4.  Preliminary  Considerations 

Using  the  mathematical  model  described  in  Section  I,  we  would  like  to  decide  whether  Ao  or  A! 
is  the  true  intensity  of  the  underlying  Poisson  process  that  gives  rise  to  our  observations  {Z(x)}. 
Clearly,  a  likelihood  ratio  test  (LRT)  is  called  for  [7,  p.  11].  Unfortunately,  a  formula  for  the  LR  of 
the  continuum  of  observations  {Z(x)}  does  not  seem  to  be  available  in  the  literature.  In  fact,  even 
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if  we  based  our  decision  in  a  finite  number  of  samples  of  {Z(x)},  say  Z(xx ),...,  Z(x^),  the  LR 
would  be  obtained  by  inverting  the  A'-variate  characteristic  function  of  Z(xt Z(x^)  under 
each  hypothesis  to  obtain  the  A'-variate  density  under  each  hypothesis.  The  quotient  of  these 
densities  would  yield  the  LR.  This  approach  is  impractical  unless  A’  is  very  small.  It  is  because  of 
these  considerations  that  we  introduce  the  following  suboptimal  approach. 

B.  A  Suboptimal  Detector 

As  a  suboptimal  detector  scheme,  we  propose  that  the  received  image  {Z(x)}  be  passed  through 
a  linear  filter  g  (to  be  chosen  later  in  Section  II- E)  and  then  sampled,  as  shown  in  Fig.  2.  The  final 
discrimination  of  hypotheses  will  be  based  on  the  sample  T.  More  precisely,  let 

F(i)  =  j g{x  -  u)Z(u)du,  (2.1) 

where  Z(-)  is  given  by  (1.1).  (Here  and  in  the  sequel  integrals  are  understood  as  being  over  all  of 
IR2  unless  otherwise  indicated.  One-dimensional  integrals  over  IR  are  indicated  by  f^).  Observe 
that  if  we  set 

</(•*•)  =  J  g(x  -  u)h{u)du,  (2.2) 

then 

Y[x)  =  £$(*-*„).  (2.3) 

Clearly  (2.3)  has  the  same  form  as  (1.1).  In  other  words,  when  the  shot-noise  process  { Z(x )} 
is  passed  through  the  linear  system  g ,  the  output  {T(i)}  is  also  a  shot-noise  process.  The  final 
processing  step  shown  in  Fig.  2  is  sampling.  We  set 

T  *  V'(0)  =  £$(-*„),  (2.4) 

V 

and  perform  an  LRT  based  on  T. 


C.  The  Likelihood  Ratio  for  T 

Let  F,(t)  =  P,(T  <  t),  i  =  0, 1,  be  the  cumulative  probability  distribution  of  T  given  that  A, 
is  the  true  intensity  of  the  underlying  Poisson  process  that  models  the  location  at  which  arriving 
photons  strike  the  imaging  system.  Clearly,  if  no  photons  arrive,  T  =  0.  The  probability  of  this 
event  is 

P,(T  =  0)  =  e"A\ 

where 

A,  =  J  A,(i )dx. 

In  our  applications  we  have  A,  <  oo.  We  thus  expect  F,(t)  to  have  the  form 


F.(t) 


[  ft(T)dr, 

J  -OO 

e_A*  +  /  /i(r)  dr, 

J  —OO 


t  <  0, 
t  >  0, 


for  some  nonnegative  /,  satisfying 

I"  f,(r)dr  =  1  -e-AV 

J  -  OO 


(2.5) 
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Remark:  In  general  F,(<)  could  have  jumps  at  points  t  ^  0,  depending  on  the  behavior  of  g.  In 
our  examples,  this  does  not  occur. 

Assuming  equally  likely  hypotheses,  if  we  observe  T  =  t,  we  perform  the  LRT 


( 

Hi 

> 

o 

1 

> 

> 

< 

0, 

t  =  0, 

Ho 

hit) 

fait) 

Hi 

> 

< 

H0 

1, 

t  ^  0. 

The  numerical  calculation  of  fo  and  f\  is  discussed  in  Appendix  A. 

D.  The  Probability  of  Error 

The  probability  of  error  incurred  using  the  LRT  (2.6)  is  denoted  by  Pe;  we  can  write  an  expression 
for  Pe  as  follows.  First,  let  Do  denote  the  set  of  t  such  that  we  decide  in  favor  of  Ho-  Clearly, 
0  6  Do  if  and  only  if  A0  <  Ai.  For  t  0,  i  6  Do  if  and  only  if  f0(t)  <  Let  D\  denote  the 

complement  of  Do.  We  write  D\  =  Dq.  Then,  under  equally  likely  hypotheses, 

Pe  =  \{<Pi[TeD0}  +  Po[T  €/?!]} 

=  +  IjdF^t)  -  dF0(t)\y  (2.7) 

Let  Ia(1)  denote  the  indicator  function  of  a  set  A.  In  other  words,  l  Ait)  =  1  if  t  £  A  and  /^(t)  =  0 
otherwise.  Using  (2.7)  and  (2.5)  we  can  then  write 

Pe  =  +  -e-*°}IDo(0)+  [  (Mt)-fo(t))dt\ 

2  l  •/£>0n{o}'  J 

=  ^{l  +  [e-A‘  -e~A°}/Do(0)  + /D'(fi(0- fo(t))dt}.  (2.8) 

In  our  applications,  Do  will  turn  out  to  be  an  interval,  or  a  union  of  disjoint  intervals.  Hence  the 
last  term  in  (2.8)  can  easily  be  computed  if  we  have  a  simple  way  to  evaluate 

Gi(t)  =  f  f,(r)dr.  (2.9) 

J  —  OO 

We  discuss  this  further  in  Appendix  B. 


E.  Selecting  the  Filter  g 


Ideally  we  would  like  to  select  g  to  minimize  the  probability  of  a 
dependence  of  Pe  on  g  is  not  readily  apparent,  we  introduce  the 
selecting  g.  We  would  like  to  choose  g  to  maximize  the  generalized 
eq.  (18)],  [2,  eq.  (6)]) 

ini  ~  Mo)2 

+  CTq 

where  /j,  and  of  are  the  mean  and  variance  of  T  under  hypothesis 
shot-noise  process,  it  follows  from  (2.4)  and  [5,  pp.  382-383]  that 


decision  error,  P..  Since  the 
following  ad-hoc  criterion  for 
“signal-to-noise  ratio”  (cf.  [1, 


(2.10) 

i  =  0,1.  Since  Y  in  (2.3)  is  a 


A 

Hi  = 


E,m 

J  g(-x)\,{x)dx. 


(2.11) 
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and 


o?  ±  E,[(T  -  m)2} 

=  J  g(-x)2Xi{x)dx. 


(2.12) 


Using  (2.11)  and  (2.12),  it  follows  from  the  Cauchy- Schwarz  inequality  that  g  maximizes  (2.10)  if 
and  only  if 

g{  ^)  ^ \  \  /  *.j  (2.13) 

Ai(x)  +  A0(x) 

where  c  is  an  arbitrary  constant.  Using  (2.2),  (2.13)  then  becomes 


J g(-x  -  u)h(u) du  = 


Ai(x)  -  A0(j) 
Ai(x)  +  A0(x) 


(2.14) 


Unfortunately,  (2.14)  may  not  have  a  solution.  For  example,  if  g  and  h  are  square  integrable,  the 
left-hand  side  of  (2.14)  will  be  a  continuous  function  of  x,  while  the  right-hand  side  need  not  be,  as 
is  the  case  in  our  examples  in  Section  III-C.  In  order  to  avoid  this  problem,  we  a  priori  constrain 
g  to  be  of  the  form 

K 

9ix)  =  ^WkHx  +  xk),  (2.15) 

k=l 

where  K  and  the  locations  xj,...,x/c  are  chosen  in  some  heuristic  fashion,  perhaps  as  a  uniform 
grid,  and  the  weights  are  chosen  to  maximize  (2.10).  Observe  that  if  (2.15)  holds,  then  it  follows 
from  (2.1)  that 

K 

Y{x)  =  5Zt»*Z(x  +  x*), 

k=\ 

and  thus 

K 

T  =  U(0)  =  5>*Z(xfc), 

k=\ 

i.e.,  the  test  statistic  T  is  a  weighted  superposition  of  the  measurements  Z(xi), . . . ,  Z(xk ).  Now,  let 
w  =  [u>i, . . .,  tu/f]'  and  Z  =  [Z(xi ),...,  Z(xK)}'.  Let  m,  =  E,[Z]  and  T,  =  E,[(Z  -  m,)(Z  -  m,)']. 
Then  T  =  w'Z,  and 

Hi  =  w'm,  and  a2  =  w'r.w,  (2.16) 

where  the  fcth  entry  of  the  vector  m,  is  (cf.  (2.11)) 


J  h{xk  ~  x)A ,(x)dx 


and  the  kl  entry  of  the  matrix  is  (cf.  (2.12)) 


j  h(xk  -  x)h(x(  -  x)A,(x)  dx. 


(2.17) 


(2.18) 


Letting  m  =  mi  -  mo  and  T  =  To  +  Tj,  we  see  that  under  the  constraint  (2.15),  (2.10)  becomes 


Iw'ml2 

wTw 


(2.19) 


By  the  Cauchy-Schwarz  inequality,  w  maximizes  (2.19)  if  and  only  if  Tw  =  cm,  where  c  is  an 
arbitrary  constant.  For  the  numerical  examples  in  Section  III-C  we  take  c  =  1,  and  Tw  =  m  is 
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easily  solved  using  the  NAG  routine  F04ASF.  The  NAG  routine  D01FCF  is  used  to  compute  (2.17) 
and  (2.18). 

Remark:  One  of  the  reviewers  has  suggested  selecting  g  so  that 

9{-x)  =  f  g(-x  -  u)h(u)du  =  ln^T^r.  (2-20) 

J  A0(x) 

the  idea  being  that  (2.4)  would  then  be  equal  (up  to  an  additive  constant)  to  the  logarithm  of  the 
LR  of  the  point  process  itself.  A  related  idea  was  used  in  [10].  Unfortunately,  (2.20)  may  not  have 
a  solution  for  the  same  reasons  given  following  (2.14). 

III.  Hypothesis  Testing  and  Error  Performance 

In  this  section  we  apply  the  preceding  ideas  to  several  examples.  For  comparison,  we  also  discuss 
the  consequences  of  assuming  that  T  has  a  Gaussian  distribution  under  each  hypothesis;  this 
assumption  was  used,  under  somewhat  higher  light  levels,  in  [4,  10]  based  on  the  Central  Limit 
Theorem.  Under  certain  conditions  this  approximation  is  adequate  [6],  and  it  avoids  the  burden 
of  computing  /;(<)  and  Gift)  by  the  numerical  evaluation  of  inverse  Fourier  transforms.  However, 
under  the  low-light-level  conditions  considered  here,  we  do  not  expect  the  Gaussian  approximation 
to  work  well,  and  this  is  indeed  the  case. 

A.  Likelihood  Ratio  Test 

Following  the  observation  T  =  t,  the  LRT  is  given  by  (2.6).  In  the  numerical  examples  discussed 
below,  we  plot  /0(<)  and  fi(t)  (Figs.  4-6)  and  see  that  the  equation 

/»(«)  =  , 

fo(t) 

has  at  most  one  solution  of  interest,  denoted  by  tj.  Hence,  if  t  ^  0,  the  LRT  reduces  to  the  single 
threshold  test 

Hi 

t  >  X], 

Ho 


B.  Gaussian  Test  (GT) 

A  simple  test  to  write  down  is  the  following.  Let 


Mt)  = 


- I _ e-(«-M.)2/2  o] 

v/2 la* 


(3.1) 


where  /i;  and  of  are  given  by  (2.16).  In  other  words,  we  are  assuming  that  T  is  normal,  but  with 
the  correct  mean  and  variance.  We  consider  the  test 

PjW  >'  , 

»>«)  H. 

In  the  examples  discussed  below,  of  >  o%,  and  the  decision  regions  for  this  test  are  easily  shown 
to  be 

D0  =  jf  €  1R  :  -7  -  -  <  t  <  7  -  -j  and  D\  =  Dq, 
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where  a  =  a2  -  Oq,  b  =  p\Oq  -  hqo\,  and 


a 

7  = 


b2 


In  the  examples  discussed  below,  the  integrals  fZ^~b^a  Pi(t)  dt  are  negligible,  and  so  we  consider 
the  following  single-threshold  test.  We  set  r?  =  7-6/a  and  use  Do  =  (-00,7?).  Even  with  this  choice 
of  Do,  the  probability  of  error  is  given  by  (2.8),  which  requires  complicated  numerical  integration 
as  discussed  in  Appendix  B.  However,  we  also  consider  Pgc  the  “Gaussian  approximation”  of  Pe, 
that  we  define  by 

Po*  =  \{f  Pii*)dt  +  J  Po(0*}- 

In  the  examples  below,  Pae  was  computed  easily  with  the  NAG  routine  S15ABF  for  evaluating  the 
cumulative  distribution  of  the  standard  normal  density. 


C.  Examples 


We  compare  the  error  performance  of  the  two  tests  in  three  examples.  In  each  example  we 
consider  three  cases,  sampling  at  A’  =  1,5,  and  10  points.  For  this  purpose,  let 

TZ  =  {x  =  (u,  v)  :  0  <  u  <  1  and  0  <  v  <  1}, 

x  x  _  /  c0Io(x)+n,  x  e  n, 

Mz)  -  {  0,  xtK, 

\  t„\  /  cl /e(x)  +  n,  x  €  U, 

hlX>  =1  0,  xfR, 

h(x)  =  e-50(u2+v2), 

where  I0(x)  and  /e(x)  are  the  indicator  functions  of  the  shaded  regions  in  Fig.  3,  and  n  is  a  constant 
background  level.  The  sampling  points,  xi,. .  .,xio,  are  shown  in  Fig.  4  and  are  explicitly  listed  in 
Table  I. 


TABLE  I 

The  Sampling  Points 


When  K  =  1,  we  set  wj  =  1  and  hence  T  =  Z(x  1). 

Example  1:  In  this  example  c0  =  ct  =  100  and  n  =  10.  Then  A0  =  43.00  and  Aj  =  44.87.  The 
weights  {u>fc}  for  the  case  K  =  5  are  1.100019, 2.603646x  10-1,  -8.424318x  10-1,  - 1.044526 x  10-2, 
-8.713632  x  10-2,  and  for  K  =  10,  1.191592,  3.038038  x  10"1,  -9.131236  x  10"\  4.417913  x  10~3, 
-7.798749 x  10“2, 1.564578 x  10-2,  -2.080950 x  10'1, 9.888259 x  10-2,  -9.350529 x  10-2,  4.023416x 
10-3.  Table  II  contains  the  means  and  variances  of  the  filtered  point  process  for  the  three  cases. 
Table  III  contains  the  thresholds  for  the  two  tests  under  consideration,  the  corresponding  decision 
regions  Do  and  probabilities  of  error.  The  value  of  Poe  is  also  included.  A  plot  of  fo  and  f\  for 
each  case  is  shown  in  Fig.  5. 


k 

Xk 

6 

(0.275,0.825) 

7 

(0.5,0.65) 

8 

(0.725,0.5) 

9 

(0.5,0.325) 

10 

(0.5,0.175) 

k 

Xk 

1 

(0.5, 0.5) 

2 

(0.7,0.175) 

3 

(0. 725,0.325) 

4 

(0.725,0.825) 

5 

(0.275,0.5) 
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TABLE  II 

Parameters  for  Example  1 


K  =  1 

cm 

O 

II 

S< 

Po 

1.52 

-1.27 

-1.38 

m 

0.422 

1.53 

1.61 

p\ 

4.40 

3.25 

3.21 

m 

2.58 

2.98 

2.99 

TABLE  III 

Test  Performance  in  Example  1 


o 

Test 

9 

Do 

Pe 

P Ge 

1 

LRT 

2.57 

0,f?) 

0.0936 

- 

GT 

2.65 

0,1?) 

0.0949 

0.0887 

5 

wmm 

0.720 

(-<*>, p) 

0.0523 

- 

GT 

0.773 

(-<*>, T)) 

0.0525 

0.0628 

10 

LRT 

0.650 

0.0506 

- 

GT 

0.706 

(-<».»?) 

0.0507 

0.0617 

From  Table  III,  we  see  that  a  reduction  of  44.1%  in  Pe  is  obtained  for  the  LRT  when  using  K  =  5 
instead  of  K  —  1,  and  a  further  reduction  of  3.25%  is  obtained  by  using  K  —  10. 

Example  2:  In  this  example  co  =  Ci  =  50  and  n  =  5.  Then  Ao  =  21.5  and  Ai  =  22.4.  The  weights 
{w*}  for  the  case  K  =  5  are  1.100019,  2.603646  x  10~\  -8.424318  x  10'1,  -1.044526  x  10~2, 
-8.713632  x  10-2,  and  for  K  =  10,  1.191592,  3.038038  x  10"1,  -9.131236  X  10"1,  4.417913  x  10~3, 
-7.798749 xl0~2,  1.564578 x  10"2,  -2.080950x  10'\ 9.888259x  10~2,  -9.350529 x  10~2,  4.023416 x 
10~3.  Table  IV  contains  the  means  and  variances  of  the  filtered  point  process  for  the  three  cases. 
Table  V  contains  the  thresholds  for  the  two  tests  under  consideration,  the  corresponding  decision 
regions  Dq  and  probabilities  of  error.  The  value  of  Poe  is  also  included.  A  plot  of  fo  and  /j  for 
each  case  is  shown  in  Fig.  6. 


TABLE  IV 

Parameters  for  Example  2 


K  =  1 

K  =  5 

O 

II 

m0 

0.759 

-0.63 

-0.691 

al 

0.211 

0.760 

0.800 

m, 

2.20 

1.62 

1.601 

A. 

1.29 

1.49 

1.50 

TABLE  V 

Test  Performance  in  Example  2 


Test 

V 

Do 

P' 

PG e 

LRT 

1.29 

Hxn 

0.177 

- 

GT 

1.44 

mmm 

0.182 

0.157 

5 

LRT 

0.389 

(-00, r?) 

0.129 

- 

GT 

0.464 

(-00,9) 

0.129 

0.138 

LRT 

0.360 

(-00,7?) 

0.126 

- 

GT 

0.423 

(-00,7?) 

0.127 

0.137 
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From  Table  V,  we  see  that  a  reduction  of  27.1%  in  Pe  is  obtained  for  the  LRT  when  using  K  =  5 
instead  of  K  =  1,  and  a  further  reduction  of  2.32%  is  obtained  by  using  K  =  10. 

Example  3:  In  this  example  Co  =  c\  =  16  and  n  =  5.  Then  Ao  =  10.28  and  Ai  =  10.58. 
The  weights  {tn*}  for  the  case  K  =  5  are  7.590372  x  10— 1 ,  1.791687  x  10— 1 ,  -5.971229  x  10~ 1 , 
-8.206410X  10"3,  -5.-151747 x  10-2,  and  for  I<  =  10,  8.767560x10-*,  2.073143x  10"1,  -6.402233  x 
10_1,  7. 654726  x  10“3,  -4.801134  x  10"2,  1.534052  x  10"2,  -1.978471  x  10"1,  7.396542  x  10"2, 
- 1.243330  x  10— 1 ,  2.736930  x  10~2.  Table  VI  contains  the  means  and  variances  of  the  filtered  point 
process  for  the  three  cases.  Table  VII  contains  the  thresholds  for  ihe  two  tests  under  consideration, 
the  corresponding  decision  regions  Do  and  probabilities  of  error.  The  value  of  P<je  is  also  included. 
A  plot  of  fo  and  f\  for  each  case  is  shown  in  Fig.  7. 


TABLE  VI 

Parameters  for  Exmaple  3 


(231 

K  =  5 

(JBLJI 

m0 

0.457 

-0.0894 

-0.122 

mm 

0.174 

0.199 

0.204 

0.918 

0.414 

0.398 

mm 

0.518 

0.305 

0.316 

TABLE  VII 

Test  Performance  in  Example  3 


B9 

Test 

7? 

Do 

Pe 

1 

LRT 

0.675 

HZEfl 

0.335 

- 

[0,7?) 

0.348 

0  310 

5 

LRT 

0.188 

(-00,7?) 

0.312 

- 

(-00,  77) 

0.313 

0.303 

10 

LRT 

0.172 

(-00,77) 

0.309 

- 

GT 

0.212 

(-00,  77) 

0.311 

0.300 

From  Table  VII,  we  see  that  a  reduction  of  6.87%  in  Pe  is  obtained  for  the  LRT  when  using  K  =  5 
instead  of  K  =  1,  and  a  further  reduction  of  0.962%  is  obtained  by  using  v  =  10. 

IV.  Discussion  and  Conclusion 

We  have  treated  image  detection  at  low  light  levels  as  a  binary  hypothesis-testing  problem  based 
on  a  1-dimensional  test  statistic.  This  statistic  was  obtained  by  filtering  the  received  image  and  then 
sampling  at  one  point.  The  filter  we  used  was  obtained  by  maximizing  an  ad-hoc  signal- to- noise 
ratio.  In  the  examples  we  considered,  we  found  that  the  largest  weights  of  the  filter  (for  the  cases 
A'  -  5  and  A'  =  10)  correspond  to  the  locations  where  the  “o”  and  the  V  do  not  overlap.  This 
makes  intuitive  sense:  at  xt  the  “e”  is  present  and  w\  is  the  largest  positive  weight;  at  13  the  “o”  is 
present  and  W3  is  the  largest-magnitude  negative  weight.  Since  ij  and  X3  are  the  most  important 
points  for  discrimination  between  H 0  and  H\,  we  observed  little  improvement  in  peformance  when 
using  A'  =  10  instead  of  I\  =  5.  The  lower  the  intensity  of  the  point  process,  the  harder  it 
is  to  discriminate  between  the  hypotheses.  The  samples  in  Example  3  bear  less  information  for 
discrimination  than  those  in  the  other  two  examples.  This  resulted  in  a  much  smaller  improvement 
in  performance  in  Example  3  than  in  the  other  two  examples  when  going  from  A  =  1  to  K  =  5.  We 
compared  our  LRT  with  a  test  that  uses  Gaussian  densities.  From  the  results  of  the  examples,  it  is 
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clear  that  the  distribution  function  of  T  is  not  Gaussian  under  either  hypothesis  (see  Figures  5-7). 
It  is  interesting  to  observe,  though,  that  the  probability  of  error  Pe  is  very  similar  for  the  two  tests, 
i.e.,  Pe  is  not  very  sensitive  to  the  value  of  77.  We  note  that  the  Gaussian  approximation  Pce  of  the 
probability  of  error  is  neither  an  upper  nor  a  lower  bound  for  Pe  since  sometimes  it  over  predicts 
(by  21.9%  in  Example  1,  K  =  10)  and  sometimes  it  under  predicts  (by  11.3%  in  Example  2,  K  =  1) 
the  value  of  Pe.  Hence,  Pee  is  not  a  reliable  quantity  for  estimating  Pe. 

Further  research  should  be  devoted  in  studying  the  behavior  of  Pe  as  a  function  of  the  weights 
{infc},  as  a  function  of  the  sampling  points  {xjt},  and  as  a  function  of  the  total  number  of  weights 
of  the  filter  g.  It  is  very  important  that  fast  and  accurate  methods  be  found  in  order  to  compute 
the  functions  /,(<)  and  G,(t),  since  straightforward  applications  of  numerical  integration  are  rather 
time  consuming. 


Appendix  A 

Evaluation  of  the  Likelihood  Ratio 


In  order  to  perform  the  test  (2.6),  we  need  to  compute  the  “density”  function  of  T  under  each 
hypothesis.  Loosely  speaking,  this  can  be  accomplished  by  computing  the  inverse  Fourier  trans¬ 
form  of  the  characteristic  function  of  T.  Several  methods  have  been  proposed  (see  [3]  and  the 
references  therein)  for  carrying  out  this  computation.  Our  numerical  solution  relies  on  the  use  of 
the  quadrature  subroutines  D01FCF  and  D01AMF  from  the  NAG  library. 

For  the  purpose  of  computing  the  LR  function  and  the  probability  of  error  for  the  test  (2.6),  we 
introduce  the  moment  generating  function  of  T,  denoted  by  Af,(s);  it  is  given  by  [5,  p.  381] 


M,(s)  =  E,(esTl 


From  (2.5)  it  follows  that 

=  exp  j 

f  A,(x)[e^-*>  -  l]dxj. 

(A.l) 

where 

M,(s)  = 

e-Al  +  Ki(s), 

(A.2) 

With  s  =  a  +  ju ,  let 

A \(s)  = 

/  e3tft{t)dt. 

J  —  OO 

(A.3) 

4  Re{j 
■  /A'(: 

(  A,(x)es®(-r)dxJ 

cos  (ug(-x))  dx, 

and 


S.^u;)  =  Im  jy 

=  J  A,(x)e<7^'x'  sin(u>(j(-x))  dx. 

We  compute  C°(ui)  and  S°(oj)  numerically  using  the  NAG  routine  D01FCF.  If  we  set  s  =  ju>  in 
(A. 3)  and  take  inverse  Fourier  transforms,  we  obtain,  since  Re{ K,(ju)eJwt}  is  an  even  function  of 


1 

fx{t)  =  -  {e~*A,-C’*w*1  cos(5,(u;)  -  ut)  -  e-A'  cos(u>t)}  do:,  (A.4) 

7T  JO 

where,  for  convenience  of  notation,  we  write  C,  and  Si  instead  of  Cf  and  Sf  when  a  =  0. 
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Appendix  B 

Evaluation  of  the  Probability  of  Error 


In  order  to  compute  Pe,  we  first  need  to  compute  (2.9).  To  this  end,  let  G,(s)  denote  the  Laplace 
transform  of  G,.  More  precisely, 


Assuming  that 


Gi(s)  =  r  Gi(t)eatdt,  Re{s)  <  0. 

J  —  OO 


lim  estGi(t)  =  0,  Re{s}  <  0, 

OO 


(B.l) 


which  is  true  in  our  applications,  integration  by  parts  yields 

Gi(s)  =  -~Kx(s). 


(B.2) 


Combining  (B.2)  with  (A. 2)  and  (A.l),  we  can  substitute  in  (B.l)  and  take  inverse  Fourier  trans¬ 
forms  (after  writing  s  =  a  +  ju,  with  a  <  0),  to  obtain 


erf.,--: f  rw*,, 

TT  J  0  (Tz  +  M* 


(B-3) 


where 


=  e  (A‘  cos( 5f  (w)  -  u;t) -f  u;sin(5,CT(u;)  -  u;<)] 

-  e~x,(cr  cos(ut)  -  u;sin(w<)). 

The  integrals  (A. 4)  and  (B.3)  are  computed  numerically  using  the  NAG  routine  D01AMF. 
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Wavelet  Transforms  for  Discrete-Time  Periodic  Signals 

John  A.  Gubner,  Member,  IEEE,  and  Wei-Bin  Chang 


Abstract —  A  tutorial  development  of  wavelet  transforms  for 
discrete-time  periodic  signals  is  presented  using  only  basic  con¬ 
cepts  from  linear  algebra.  Like  the  discrete  Fourier  transform, 
these  wavelet  transforms  are  shown  to  be  unitary  operators  that 
exhibit  frequency-localization  properties.  They  are  also  shown 
to  exhibit  time-locaii2ation  properties  and  to  be  well  suited  for 
signal-compression  applications.  Wavelet  transforms  for  im¬ 
ages.  i.e..  matrices,  are  briefly  discussed.  Wavelet  transforms 
are  constructed  that  are  fast  in  the  sense  that  their  complex¬ 
ity  is  no  greater  than  that  of  the  corresponding  fast  Fourier 
transform. 

I.  Introduction 

The  study  of  wavelets  has  drawn  considerable  attention 
from  the  signal  processing  community  since  the  appearance 
of  the  papers  by  Daubechies  [2,  3]  and  Mallat  [8,  9,  10]. 
I  nfortunately,  a  complete  understanding  of  their  work  re¬ 
quires  a  solid  background  in  functional  analysis.  However, 
just  as  the  discrete  Fourier  transform  for  finite-length  data 
records  can  be  rigorously  derived  using  only  concepts  from 
linear  algebra,  we  do  the  same  for  wavelet  transforms  in 
this  tutorial.  It  is  hoped  that  by  restricting  attention  to 
this  finite-dimensional  framework,  the  presentation  will  be 
more  accessible. 

There  are  certain  differences  between  continuous-time, 
infinite-dimensional  wavelets  and  the  finite-dimensional, 
discrete-time  wavelets  considered  here.  For  example,  the 
analog  of  =  \/2  t/'j +  i.o(2<)  [2,  p.  958]  does  not  hold 

here.  Another  difference  is  that  in  the  infinite-dimensional 
case,  one  deals  with  a  single  pair  of  operators,  G  and  H , 
while  for  discrete-time  periodic  signals  we  have  a  sequence 
of  operators,  (7;  and  .  Furthermore,  instead  of  expres¬ 
sions  such  as  (G*)nr  [2,  p.  947],  we  consider  expressions 
of  the  form  Go  G^l^lz\  in  the  former  case,  there  is  the 
recursion.  (G*)nz  =  G*((G*)n_1i],  while  for  us  there  is  no 
simple  recursion  to  obtain  Gg  ■  GJJ_]2  from  GJ  G^_2i- 
Fortunately,  these  differences  do  more  to  reveal  the  struc¬ 
ture  of  wavelet  transforms  than  to  complicate  our  devel¬ 
opment. 

Tutorials  that  discuss  continuous-time  wavelets  and  re¬ 
lated  topics  can  be  found  in  [5,  6,  13,  15].  Extensive  bib¬ 
liographies  are  contained  in  (5,  6,  13], 


This  work  was  *upport*d  in  part  by  the  Air  Force  Office  of  Scien¬ 
tific  Research  under  Grant  AFOSR-90-0181 . 

The  authors  are  with  the  Department  of  Electrical  and  Computer 
Engineering.  University  of  Wisconsin.  Madison,  WT  53706. 


A.  What  are  Wavelets? 

Let  Xo  denote  the  space  of  discrete-time  signals  xo (n) 
with  period  N,  where  N  is  a  power  of  2.  Before  defining 
wavelets,  recall  that  from  the  theory  of  discrete  Fourier 
transforms,  every  signal  xo  €  Xo  can  be  expressed  as 

N-l 

*o(n)  =  ^2(x0,  Ek)Ek(n), 

k~0 

where 

£*(")  =  3»eJ*kn’ 

N-l 

(■Eo.yo)  =  *°(R)y°(n), 

n=0 

and  the  overbar  denotes  the  complex  conjugate.  The  rea¬ 
son  the  decomposition  works  is  that  the  signals 
form  an  orthonormal  basis  for  Vo.  While  direct  computa¬ 
tion  of  the  Fourier  coefficients, 

1  N~‘ 

(x°,£i)  =  -j=  £  zo(n)e-]#kn, 

for  k  =  0, . . . ,  N  —  1  would  require  JV2  multiplications,  the 
fast  Fourier  transform  can  compute  these  coefficients  with 
order  N  log2  N  multiplications. 

In  this  paper  we  consider  orthonormal  “wavelet”  bases 
for  A’o  consisting  of  signals  and  For  vJt, 

j  =  1 . J  and  for  each  j,  i  -  0,...,N/23  -  1.  For 

ffiji,  i  =  0,...,N/2J  —  1.  Since  these  signals  form  an 
orthonormal  basis,  every  xo  €  A’o  can  be  written  in  the 
form 

J  ,N/V-l 

*o(")  =  51  (  Ye  <*o,  V'ii)V'ji(n) ) 

y=l  .=0 

N/2J  —  l 

+  ^2  i*o 

»=o 

The  signals  are  called  wavelets,  and  the  signals  w,  are 
called  scaling  functions;  examples  for  the  case  N  =  128 
and  J  =  5  are  shown  in  Figs.  1-6.  The  numbers  (xo ,4’jt) 
and  {xo,ifiji)  are  called  wavelet  coefficients,  and  the  fast 
wavelet  transform  is  an  algorithm  for  computing  them  with 
fewer  than  N1  multiplications;  in  some  cases  this  can  be 
done  with  order  JV  log2  N  multiplications. 

We  now  compare  and  contrast  the  sinusoidal  basis  func¬ 
tions  Ek  with  the  wavelets  Wi  and  pj, .  The  constant  func¬ 
tion,  Eo(n)  =  l/VN,  exhibits  no  variation  at  all,  while  the 


HO 


signals  exhibit  the  least  variation  of  all  the  wavelet 

signals.  At  the  other  extreme,  the  signal  En- i(n)  exhibits 
the  most  variation  of  all  the  sinusoids  £7*(  v»),  while  the  sig¬ 
nals  t/’i  ,(n)  exhibit  the  most  variation  of  all  the  wavelet 
signals.  One  important  difference  between  the  wavelets  in 
Figs.  1-6  and  the  sinusoids  Ek  is  that  the  are  zero 

for  most  values  of  n.  For  example,  in  Fig.  1,  we  see  that 
tin .2<n)  is  nonzero  only  for  n  =  2, 3, 4, 5.  Thus,  looking 
at  wavelet  coefficients  can  help  us  localize  high-variation 
features  in  time  rather  precisely.  Slower-varying  signals 
cannot  be  localized  as  well. 

B.  Outline  of  the  Paper 

The  remainder  of  the  paper  is  organized  as  follows.  In 
Section  (I  we  focus  on  the  basic  properties  of  pyramid- 
type  algorithms  [1],  which  will  be  the  key  to  defining  and 
computing  wavelet  transforms  in  Section  III.  In  particu¬ 
lar.  we  show  that  such  algorithms  are  natural  vehicles  for 
frequency  localization  and  for  signal  compression.  We  also 
show  how  these  algorithms  interact  with  shift  operators  to 
exhibit  time-localization  behavior. 

In  Section  III  we  define  wavelet  transforms  for  discrete¬ 
time  periodic  signals.  We  then  define  the  corresponding 
wavelet  bases  in  Section  III- A.  In  Section  III-B  we  dis¬ 
cuss  fast  wavelet  transforms  whose  complexity  is  no  greater 
than  that  of  the  corresponding  fast  Fourier  transform.  Nu¬ 
merical  examples  are  presented  in  Section  III-C,  where 
we  discuss  frequency  localization,  time-frequency  filtering, 
and  signal  compression  using  discrete  wavelet  transforms. 
In  Section  III-D  we  briefly  sketch  an  extension  of  wavelet 
transforms  to  two-dimensional  signals  such  as  images. 

In  Section  IV  we  extract  the  key  ideas  from  [2,  Sec¬ 
tions  3.  A.  and  4.B.]  to  give  a  self-contained  derivation  of 
the  sequence  of  operators  used  to  implement  fast  wavelet 
transforms  for  discrete-time  periodic  signals. 

II.  Pyramid  Algorithms 

The  basic  building  block  for  pyramid-type  algorithms  is 
a  mapping  from  an  input  space  X  into  a  product  space, 
say  V  x  V ,  as  discussed  in  Section  II- A .  In  Section  II-B 
we  iterate  a  sequence  of  these  operators  to  obtain  a  De¬ 
composition  Algorithm  and  a  Reconstruction  Algorithm. 
Applications  to  signal  compression  and  time  localization 
are  discussed  in  Sections  II-C  and  II- D,  respectively. 


Now,  if  X  is  a  third  inner-product  space,  and  if  G.  X  — ■  U 
and  H:  X  — ►  V ,  are  linear  operators,  we  define  the  map¬ 
ping  T.X  —  U  x  V  by 


Next,  if  the  adjoint  operators  G’ :  U  — *  X  and  Hm :  V  —  ,Y 
exist,  then  the  adjoint  T’:U  x  V  —*  X  also  exists  and  is 
given  by 


[G'  //'] 


G'u  +  H‘v. 


Our  interest  is  in  isometries;  i.e. ,  operators  that  preserve 
inner  products,  or  equivalently,  operators  for  which  T'T  = 
/.  For  such  operators,  x  can  obviously  be  recovered  from 
Tx. 

Remark:  If  T'T  =  /  and  T  is  onto,  then  TT'  =  1  as 
well  [4,  p.  185].  Such  an  operator  is  said  to  be  unitary.  Of 
course,  if  dim  A'  =  dim  U  x  V  <  oo,  then  T'T  =  I  always 
implies  TT'  =  /,  since  in  this  case  T  is  1-to-l  if  and  only 
if  T  is  onto  [7,  p.  81], 

For  future  reference  note  that  since  T  has  the  form  in 
(2),  T'T  =  /  is  equivalent  to 


G'G  +  HmH  =  /,  (3) 

and  TT *  =  /  is  equivalent  to  the  following  three  equations, 


GG*  =  I 
HH'  =  I 
GH '  =  0. 

To  summarize,  if  we  “decompose”  x  with 

u  =  Gx, 
v  =  Hx, 


(4) 

(5) 

(6) 


(7) 


and,  if  (3)  holds,  then  we  can  reconstruct  x  from  u  and  v 
with 

x  =  G'u  +  H'v. 

Remark:  As  shown  in  Section  IV-C,  in  the  wavelet  set¬ 
ting,  G  can  be  viewed  as  a  low-pass  filter,  and  H  can  be 
viewed  as  a  high-pass  filter.  Thus,  in  (7),  we  think  of  t;  as 
carrying  the  high-frequency  detail  in  x,  and  we  think  of  u 
as  carrying  the  smoother,  low-frequency  information  in  x. 


.4.  Isometries  and  Product  Spaces 


Let  (’  and  V  be  inner-product  spaces  with  respective  in¬ 
ner  products  (  ,  )  and  (•,•).  This  induces  an  inner  product 
on  V  x  T  by  taking 


([:]•[:]) 


=  (u,u)  +  (e,t>). 


B.  The  Genera/  Algorithms 

Let  Xo,X\,. .  .,X  j  and  V\ , . .  .,Vj  be  two  sequences  of 
inner-product  spaces.  Denote  the  inner  product  on  A'; 
by  (-,  )j  and  the  inner  product  on  Vj  by  For  j  = 

0 . J  -  1,  let  Gj  .Xj  —  Xj+i  and  let  Hj  .  Xj  —  Vj+i. 

Consider  the  following  procedure. 


Decomposition  Algorithm 

Let  xo  €  Xo  be  given, 
for  j  =  0  to  J  —  1 

*j+i  =  GjXj 

vj+ 1  —  Hjxj 
next  j 
end 

If  (3)  holds  for  each  pair  Gj  and  Hj ,  then  x0  can  be 
recovered  from  V\ , . . . .  vj ,  xj  by  the  following  procedure. 

Reconstruction  Algorithm 

Let  i'i , .  .  . ,  vj,  xj  be  given, 
for  ./  =  J  -  1  to  0 

x.  -  G‘x}  +  \  +  HJ  v]  + 1 
next  j 
end 

If  we  let 

V  =  V'i  x  •  x  Vj  x  Xj, 

and  if  i’i . vj,  and  xj  are  generated  by  the  Decomposi¬ 

tion  Algorithm,  then  setting 

Dxn  =  (v\ . vj,xj)  (8) 

defines  a  mapping  D  Xo  — *  Y.  We  equip  Y  with  the 
induced  inner  product,  denoted  by  (  ,  ■).  Note  that  we  are 
now  using  row-vector  notation  for  product-space  vectors. 

Remark:  We  can  extend  the  ideas  in  the  remark  at 
the  end  of  Section  II- A  to  D  by  saying  that  tq  contains 
the  highest  frequency  information  in  xo,  while  vj  and  xj 
contain  the  lowest  frequency  information  about  xo- 

Proposition  1:  Assuming  only  that  G*  and  HJ  exist,  the 
vector  xn  generated  by  applying  the  Reconstruction  Algo¬ 
rithm  to  an  arbitrary  vector  y  =  (t?i , . . . ,  Vj,  Xj)  £  Y,  is 
equal  to  D'y\  i.e.,  the  Reconstruction  Algorithm  imple¬ 
ments  D*. 

Proof:  Let  xo  be  any  element  of  Xo,  and  let  y 

be  any  element  of  Y.  We  must  show  that  if  xo  is  ob¬ 
tained  by  applying  the  Reconstruction  Algorithm  to  y, 
then  (Dx0.y)  =  (x0,x0)o-  Write 

j 

(Dx0,y)  =  +  ( xj,xj)j ,  (9) 

j=t 

where  i 'j,  .  ,!•/,  and  ij  are  obtained  by  applying  the 
Decomposition  Algorithm  to  xo-  In  addition  to  xo,  let 
i|,  ix-i  also  be  obtained  by  applying  the  Reconstruc¬ 
tion  Algorithm  to  y  Then  for  each  j  =  J, ... ,  1,  substitute 


the  following  identity  into  (9): 

+  (xhxi)i 

=  +  {Gj.  xxj.ux,), 

=  (xj  - 1 1  xj  - 1 );  - 1  • 

□ 

Corollary  2:  If  (3)  holds  for  each  pair  Gj  and  H},  then 
D’D  =  /;  i.e.,  D  is  an  isometry. 

Proof:  If  we  had  assumed  that  y  =  Dx o  in  the  pre¬ 
ceding  proof,  then  we  would  have  obtained  (Dxq,  Dio)  = 
{x0,xo)ol  i  e.,  D  preserves  inner  products.  □ 

Corollary  3:  If  (4)-(6)  hold  for  each  pair  Gj  and  H}, 
then  DD *  =  I.  If  (3)  also  holds,  then  D  is  unitary. 

Proof:  Fix  a  y  =  (vi, . . . ,  vj ,  xj)  €  Y  and  let 

xj_i,...,x o  be  generated  by  the  Reconstruction  Algo¬ 
rithm.  Now  apply  the  Decomposition  Algorithm  to  x0 
Since  x0  =  GqX i  +  HqV i  ,  Gox0  =  x\  by  (4)  and  (6).  Simi¬ 
larly,  for  j  =  1, . . .,  J  -  1,  Gj  Xj  =  Xj  +  i .  In  addition,  since 
GjHJ  =0  implies  H}G’  =0, 

Hjij  =  Hj(G’xj+i  +  HJvj+i)  =  vj+ 1. 

Hence,  DD'y  =  y.  □ 

Proposition  4:  If  D  is  unitary,  then  we  can  write  A"o  as 

Xo  =  Wx  ©  •  ■  -©  Wj  ®  Bj, 

where 

W;  =  {D*y  :y=( 0 . 0,  Vj,0, . .  .,0),Uj  G  V3) 

and 

Bj  =  {D'y  :  y  =  (0, . . . ,  0,  xj),  xj  €  Xj  } 

are  orthogonal  subspaces  of  Xo- 

Proof:  Since  D‘ D  =  /,  we  can  write  xo  =  D'{Dx o). 
and  by  (8)  it  is  clear  that  Xo  is  a  sum  of  the  indicated 
subspaces.  Since  DD *  =  /,  it  is  easy  to  see  that  Xo  is  a 
direct  sum  and  that  the  subspaces  are  orthogonal.  □ 

C.  Application  to  Signal  Compression 

Given  io  €  Xo,  apply  the  Decomposition  Algorithm  to 
obtain  Dx o  €  Y  as  in  (8).  Assuming  that  (3)  holds  for 
each  Gj  and  Hj ,  D  preserves  norms,  and  thus 

Ikollo  =  l!£>*ollJ 

=  ElMj  +  IMlS- 

;  =  i 


To  store  a  compressed  version  of  xo,  we  save  only  those  vj 
and  xj  of  significant  energy  relative  to  ||*o|lo-  To  recover 
an  estimate  of  xq ,  we  apply  the  Reconstruction  Algorithm, 
replacing  the  omitted  data  with  zeros.  Numerical  examples 
are  discussed  in  Section  III-C. 

D.  Interaction  with  “Shift”  Operators — Time  Localization 

Fix  an  integer  k,  1  <  k  <  J ,  and  consider  a  sequence  of 
operators  So,  •  .  . ,  St  such  that  Sy :  Ay  — »  Xj  and  such  that 

G,S j  =  S]  +  ,Gj,  j  =  0 . k-  1.  (10) 

Remark:  When  this  setup  is  applied  to  wavelets,  we 
take  Xj  to  be  the  set  of  column  vectors  of  length  jV/2;, 
where  ,V  is  a  power  of  2.  We  interpret  So  as  a  circular 
shift  of  2*  units,  and  Sy  as  a  shift  of  2k~J  units.  Thus  Sk 
would  be  a  shift  of  2°  =  1  unit.  See  Section  III  for  further 
discussion. 

Fix  Xo  €  A'o,  and  let  xj.t'i,  ....  xj,vj  be  given  by  the 
Decomposition  Algorithm.  Set  x0  =  SoXo  and  let  xi.Ci, 
ij  l'j  he  generated  by  applying  the  Decomposition 
Algorithm  to  xo-  Then  it  follows  by  induction  that 

*1=5}*,,  J  =  0 . k.  (11) 

At  this  point  we  restrict  our  setup  and  require  that  the 
space  + ,  =  ,YJ  +  [,  j  =  0. . . . ,  7  -  1,  so  that  S;  +  i  can 
operate  on  elements  of  Vy  +  i  =  AJ+i.  We  assume  that  in 


addition  to  ( 10), 

HjSj  =  Sj  +  iHj,  j  =  0 . k-  1.  (12) 

With  ij ,  ij  and  x; ,  £>y  as  before,  we  have,  for  j  =  0 . k  — 

1 , 

l’y+i  =  H  jij 

-  HjSjXj,  by  (11), 

=  +  by  (12), 

—  •5;  +  tt’y  +  i-  1 13) 


Now  consider  an  x0  such  that  applying  the  Decomposi¬ 
tion  Algorithm  to  Xo  yields 

Dx o  =  (0, . . .  ,0,  vt,0, . . .  ,0). 

We  claim  that  if  (3)  holds  for  each  pair  G;  and  Hj,  then 

Dx0  =  D(SoXq)  =  (0, . .  .,0,5*»»,0, ..  .,0). 

By  (13)  i)i . i  are  all  zero  and  C*  =  S*v*.  If  it  =  J, 

simply  note  that  ij  =  0  by  (11).  If  k  <  7,  the  Reconstruc¬ 
tion  Algorithm  shows  (on  account  of  (3))  that  xj, . . .  ,xk 
are  all  zero.  Then  by  (11),  x*  =  0.  Now,  by  the  Decom¬ 
position  Algorithm.  x*+i .  i>+1,  *j,  vj  are  all  zero. 

In  general  we  have 

D(S‘0x  0)  =  (0 . 0,  S[vk,0, . . 0). 

Remark:  If  D  is  unitary,  we  see  that  So  maps  Wk  into 
itself  If  k  =  7,  So  also  maps  A'j  into  itself. 


III.  Discrete  Wavelet  Transforms 

To  make  the  foregoing  discussion  concrete,  we  apply  it 
to  the  following  situation.  Let  N  be  a  power  of  2,  and  set 

Xj  ~  CN/ 2’,  7  =  0 . 7, 

where  2  <  N/2J  and  CN/2’  is  equipped  with  the  usual 

Euclidean  inner  product.  We  take  V}  =  Ay,j  =  1 . 7, 

as  was  done  in  Section  II-D. 

Convention:  As  in  the  theory  of  discrete  Fourier  trans¬ 
forms,  we  do  not  distinguish  between  vectors  in  CN/7’  and 
their  doubly-infinite  iV/2^ -periodic  extensions.  Recall  that 
the  Af-length  discrete  Fourier  transform  (DFT)  of  a  se¬ 
quence  x  €  CN  is  defined  by 

N-l 

n =0 

for  all  integers  u. 

Definition  5  (Discrete  Wavelet  Transform ):  If  the  oper¬ 
ators  Gj  .  Xj  —  Ay  +  i  and  H } :  *y  -  V;  +  1  have  the  form 

N/  2»-l 

(G;x)(m)  =  J2  9j(n  ~  2m)x(n),  (14) 

n=0 
N/23  —  1 

(Hyx)(m)  =  5Z  hy(n  -2m)x(n),  (15) 

n=0 

where  gj,  hj,  and  x  are  Af/2-7-length  periodic  sequences, 
and  where  gj  and  hj  are  such  that  (3)— (6)  hold  for  G ,  and 
Hj ,  then  Dx o  defined  by  using  these  operators  in  the  De¬ 
composition  Algorithm  is  called  a  discrete  wavelet  trans¬ 
form  (DWT)  of  *0. 

For  operators  of  the  form  (14)  and  (15),  if  1  <  fc  <  7, 
and  if 

(Syx)(n)  =  x(n  —  2*_,))  j  =  0, 

where  x  €  Ay  =  CN^’ ,  then  the  periodicity  convention 
can  easily  be  used  to  show  that  (10)  and  (12)  hold. 

Remarks:  ( i )  From  the  derivation  in  Section  1V-A  and 
from  (57)  in  the  Appendix,  it  can  be  seen  that  if  G}  and 
Hj  satisfy  (3)-(6),  then  gj  and  hj  are  discrete-time  pe¬ 
riodic  analogs  of  conjugate  quadrature  filters  in  [14] .  In 
Section  IV  we  construct  py  and  hj  under  two  additional 
constraints.  The  first  is  that  jy(n)  and  hy(n)  be  zero  for 
most  values  of  n.  The  second  is  that  Daubechies’  con¬ 
straint  (51)  hold. 

(ii)  The  Decomposition  and  Reconstruction  Algorithms 
can  be  implemented  by  analysis  and  synthesis  filter  banks, 
respectively. 

(m)  The  operators  Gy  and  Hj  in  (14)  and  (15)  can  be 
identified  with  (N/2*+l)x(N/V)  matrices;  for  an  example 
with  j  =  0,  see  (28)— (29)  below. 


i)  3 


.•4  Wavelet  Bases 


B.  Fast  Wavelet  Transforms 


We  now  construct  a  special  orthonormal  basis  for  A'o. 
Recalling  Proposition  4,  it  suffices  to  construct  orthonor¬ 
mal  bases  for  W}  and  Bj. 

For  each  space  ClV^2  ,  let  6j  denote  the  iV/2^-periodic 
Kronecker  delta,  and  set  6Jt (n)  =  6;(n  —  «).  Then  the 
Sj,.  i  =  0. . . . ,  .V/2J  -  1,  constitute  the  standard  basis  for 
=  A";  =  Vj .  Next,  set 

d},  =  (0 . 6ji . 0),  j  =  1 . J. 

where  is  in  the  jth  position  out  of  J  +  1  positions  Set 

eji  =  (0 . 0. 6ji). 

If  we  let 

tji  —  D  dj,  and  •PJi  D  ej 
and  if  we  assume  that  D  is  unitary,  then 

{<-'j0 . } 

is  an  orthonormal  basis  for  W} ,  and 

{rJ 0 . -fj.X/21-  l  } 

is  an  ort  honormal  basis  for  Bj .  The  union  of  all  these  bases 
is  called  a  wavelet  basis  for  A'o  The  signals  i’},  are  called 
wavelets,  and  the  signals  ^ j ,  are  called  scaling  functions. 
Kxamples  are  discussed  in  Section  III-C. 

Remarks:  (i)  Recalling  the  expansion  (1),  observe  that 
(■ro.  <i’j>)o  =  (ra.D‘d}l)0 

=  ( Dio.dj,) 

=  t',(«) 

Similarly,  (•r0,-pvl)n  =  x j(i) 

( ii )  By  the  results  of  Section  II-D,  we  see  that  for  k  = 

1.  ...J. 


Let  p  be  an  even  integer  such  that  2  <  p  <  .V /2J . 
As  will  be  shown  in  Section  IV,  one  can  find  a  sequence 
9(0), ....  g(p  —  1 )  such  that  if 


and  if 


0, 


n  =  0 . p-  1, 

n  =  p . V/2;  -  1, 


(19) 


M")  =  (-1)"«(1  -  n).  "=0 . JV/2>-l.  (20) 


where  j  =  0 . J ,  then  the  operators  G}  and  given 

by  (14)  and  (15),  respectively,  satisfy  (3)-(6).  The  discrete 
wavelet  transform  is  then  implemented  by  the  Decomposi¬ 
tion  Algorithm,  and  the  inverse  transform  is  implemented 
by  the  Reconstruction  Algorithm. 

We  now  show  that  for  2p  <  log2  ;V,  such  a  discrete 
wavelet  transform  is  no  more  difficult  to  compute  than 
an  A-point  fast  Fourier  transform,  which  requires  order 
.V  log 2  A  complex  multiplications.  Since  g3(n)  and  h}(n) 
have  only  p  nonzero  terms,  the  number  of  multiplications 
required  to  compute  ( 14)  and  ( 15)  for  m  =  0, . . . ,  iV/2J  +  1  — 
1  is  pN/23 .  To  implement  the  Decomposition  Algorithm 
requires  that  we  do  this  for  j  —  0, . . . ,  J  —  1 .  Hence,  the 
total  number  of  multiplications  is  always  less  than  2 pN . 
If  2 p  <  log2  V,  then  2 pN  <  N  log2  N.  Note  that  even  if 
2 p  >  log2  N ,  we  can  still  use  the  inequality  p  <  N/ 2J  to 
bound  the  number  of  multiplications  by  N7/ 2J~1,  which 
is  less  than  the  N 2  multiplications  required  to  compute 
directly  all  of  the  A-length  inner  products  appearing  in 
(1). 

In  fact,  further  speed-ups  can  be  achieved  by  adapting 
some  of  the  ideas  in  Rioul  and  Duhamel  [12,  Section  III]. 
For  example,  at  each  interation  j,  the  sums  in  (14)  and 
(15)  can  be  split  into  sums  over  n  =  even  and  n  =  odd 
This  results  in  the  sum  of  two  convolutions.  Depending  on 
the  relative  values  of  p  and  log2(jV/2^  +  1),  it  may  be  more 
efficient  to  do  this  with  fast  Fourier  transforms.  For  small 
p,  Rioul  and  Duhamel  suggest  the  use  of  short-length  fast 
running  FIR  algorithms  [11,  16]. 


Uki(n-I2k)  =  il>ki+i(n),  1  =  0,...,  N/2k  -  1,  (16) 

and 

r-j,(n-l2J)  =  spj  ^i(n),  1  =  0,...,  N/2J  -  1.  (17) 

In  particular,  it  follows  that  the  subspace  W*  (resp.  Bj) 
is  closed  under  cyclic  shifts  of  2*  (resp.  2J)  units. 

(m)  The  Reconstruction  Algorithm  shows  that  ^tl  = 
// *  A  - , .  v2,  =  and 


C.  Numerical  Examples 

To  obtain  the  g(n)  mentioned  at  the  beginning  of  the 
previous  subsection,  it  is  necessary  to  solve  (50)— (51 )  in 
Section  IV.  We  take  p  =  4  so  that  these  equations  can 
easily  be  solved  by  hand.  In  the  course  of  the  calcula¬ 
tions,  there  are  two  places  where  square  roots  must  be 
taken.  Thus,  there  are  four  possible  solutions,  depending 
on  whether  positive  or  negative  roots  are  used.  If  the  neg¬ 
ative  root  is  always  taken,  we  obtain 


9(0)  =  .482962913145 
9(1)  =  836516303738 
9(2)  =  .224143868042 
9(3)  =  —129409522551 


These  are  the  first  four  entries  in  [2,  Table  1,  p.  980]  (where 
they  are  called  h(n)).  If  we  had  always  taken  the  positive 
square  root,  we  would  have  obtained  g(n)  =  g( 3  -  n), 
n  =  0. 1.2,3.  The  two  remaining  solutions  are  -g(n)  and 
-g(n). 

Let  g  be  given  by  (21)  and  define  gj  and  h ,  by  (19)  and 
(20)  with  X  -  128  and  J  —  5.  Note  that  since  g}  and  h3 
are  real,  if  xo  is  real,  so  is  Dx q. 


most  of  the  energy  of  the  high-frequency  signal  is  in  tq . 
Keeping  in  mind  that  tq  is  a  subsampled  signal  of  length 
64  =  128/2,  and  that  16/2  =  8  and  that  63/2  truncates  to 
31,  we  set  tq(8), . . . ,  tq(31)  to  zero.  If  we  inverse  DVVT. 
we  obtain  the  signal  shown  in  Fig.  11. 

To  remove  the  low-frequency  signal  during  a  specified 
time  interval  would  be  more  difficult  because  the  time  res¬ 
olution  of  Vj  for  larger  j  is  not  as  fine. 


Example  l  (Wavelet  Basis):  The  signals  V'i2.  tf'2,2. 

<•3  2  t'4.2.  U’5,2'  anrl  *?sq  are  shown  in  Figs.  1-6,  respec¬ 
tively  The  remaining  basis  signals  can  be  obtained  by 
translation  as  indicated  in  (16)  and  (17). 

From  the  figures,  one  might  suspect  that 

f  v/2t/J  +  i.2(2n),  eJ  +  12(2n)  ^  0, 

I  U,  otherwise. 

However,  closer  inspection  reveals  that  this  cannot  be  true. 
To  see  this,  note  that  from  Figs.  1  and  2,  tr2  2  has  10 
nonzero  components,  while  cq  2  has  only  4  nonzero  com¬ 
ponents.  To  more  fully  understand  what  is  happening,  it 
is  necessary  to  consider  the  limiting  continuous-time  case 
as  in  [2]. 

Erample  2  ( Frequency  Localization):  As  a  low-frequen¬ 
cy  signal  we  consider 

Tnj,  n  =  0 .  127,  (22) 

and  as  a  high-frequency  signal  we  consider 


=  »*(  ± 


Example  4  ( Signal  Compression):  We  now  apply  the 
ideas  of  Section  II-C  to  the  DWT  shown  in  Fig.  9.  If 
we  set  tq  =  0;  i.e.,  we  set  the  first  64  elements  in  the  fig¬ 
ure  to  zero  for  a  2-to-l  compression  ratio,  and  if  we  apply 
the  Reconstruction  Algorithm,  we  obtain  the  waveform  in 
Fig.  12. 

In  Fig.  10  we  set  t>2  =  0;  i.e.,  we  set  elements  64-95  to 
zero  for  a  4-to-3  compression  ratio.  Applying  the  Recon¬ 
struction  Algorithm,  we  obtain  Fig.  13. 

As  a  final  example  of  signal  compression,  we  consider 
the  Gaussian  pulse, 


x0(n)  -  exp 


n  =  0 . 127, 


(24) 


shown  in  Fig.  14.  Its  DWT  is  shown  in  Fig.  15.  If  we  set 
both  tq  and  t/2  to  zero,  for  a  compression  ratio  of  4-to-l. 
and  apply  the  Reconstruction  Algorithm,  we  obtain  the 
result  shown  in  Fig.  16. 


We  close  this  subsection  by  noting  that  improved  results 
can  be  obtained  at  the  expense  of  using  a  larger  value  of 
P 


Xn(n)  =  cos  (  -^53n 

1  40 


n  =  0, . . . ,  127. 


(23) 


I  he  signals  are  plotted  in  Figs  7  and  8,  and  their  DWTs 
are  shown  in  Figs  9  and  10. 

Sot e:  In  the  current  setup,  if  y  =  Dx o,  then 
y  <E  IR64  x  IR32  x  IR16  x  BT*  x  IR4  x  IR4 


Hence,  when  examining  graphs  of  y(n),  one  should  keep  in 

rmnd  that  y( 0) . y( 63)  correspond  fo  tq(0),  ....  ui(63), 

y(64).  . .  y(95)  correspond  to  v2(0),  ...  t'2(31),  and  so 

on 

Clearly,  most  of  the  energy  in  the  low-frequency  signal 
is  located  in  tq  and  r4,  while  most  of  the  energy  in  the 
high-frequency  signal  is  located  in  tq . 


D.  Discrete  Wavelet  Transforms  for  Images 

We  now  briefly  sketch  a  simple  extension  of  one¬ 
dimensional  wavelet  transforms  to  handle  two-dimensional 
images;  i.e.,  where  before  xq  was  an  N  x  1  vector,  we  now 
suppose  that  zo  is  an  N  x  X  matrix,  where  N  and  X  are 
both  powers  of  2.  Let  G  and  H  be  N/ 2  x  N  matrices 
that  satisfy  (3)-(6).  Similarly,  let  G  and  H  be  X  / 2  x  X 
matrices  that  also  satisfy  (3)-(6).  Set 


T  = 

'  G  1 

and  T  = 

'  G 

H  j 

H 

By  our  assumptions  on  G,  H  and  G,  H ,  we  have  T'T  =  I, 
TT *  =  /,  T'T  =  /,  and  TT*  =  /.  Thus,  the  mapping 
xo  — >  TxqT‘  is  a  unitary  operator. 

Now  observe  that 


Example  3  ( Time- Frequency  Filtering):  Let  xo  denote 
the  sum  of  the  two  signals  in  (22)  and  (23).  We  wish 
to  remove  the  high-frequency  signal  during  the  time  inter¬ 
val  n  -  16,  .63  Let  ( tq  .....  1-5,  X5)  denote  the  DWT  of 

x.i  Considering  the  sum  of  the  DWTs  in  Figs.  9  and  10, 


Tx0f' 


Gx  qCj'  GxqH * 
Hx0G‘  Hx0H * 


(25) 


Let  xi  =  GxqG *  and  let  tq  denote  the  remaining  blocks 
of  the  matrix  in  (25).  It  should  now  be  clear  how  to  write 


(32) 


the  Decomposition  and  Reconstruction  Algorithms  using 
a  sequence  of  matrices  Gj,  Hj,  Gj,  and  Hj  that  satisfy 
(3)-(6).  If  these  matrices  are  obtained  as  in  Section  III-B, 
then  we  require  that  2  <  p,  p  <  min{W,  N}/2J . 

IV.  Construction  of  G ,  and  H,  for  Fast 
Wavelet  Transforms 

As  we  saw  in  Section  III-B,  the  key  to  fast  wavelet  trans¬ 
forms  is  to  find  g}  such  that  most  of  the  g;(n)  are  zero. 
Then  by  (20)  most  of  the  h}  (n)  will  also  be  zero.  Of  course, 
we  also  need  to  show  that  (3)-(6)  hold  for  the  correspond¬ 
ing  operators  G}  and  H}. 

Our  plan  is  to  construct  G}  from  Go  and  then  to  con¬ 
struct  Hj  from  Gj. 

A  Obtaining  H}  from  G} 

Without  loss  of  generality  we  may  work  with  Go  and  Ho- 
To  simplify  the  notation,  we  drop  the  subscript  0.  and  we 


let  M  =  .V/2.  Then  (14)  and  (15)  become 

(Gx)(m) 

.V-  1 

=  y  g(n  -  2m)x(n). 

(26) 

n  =  0 

.V-l 

(Hx)(m) 

=  y  h(n  -  2m)jr(n), 

n=  0 

(27) 

N-i 

^2  g(n  -  2m)h(n)  =  0, 

n=0 

where  the  overbar  denotes  the  complex  conjugate,  and 
6  denotes  the  Af-periodic  Kronecker  delta.  If  we  de¬ 
fine  the  A/-periodic  sequences  g<?(k)  =  g(2k),  gQ(k)  = 
g(2k  +  1  ),ht(k)  =  h(2k),  and  h0(k)  =  h(2k  +  1),  then 
(30)— (32)  become 

M-l 

y  9*{k-  m)gt(k)  +  g0{k  -  m)g0(k)  =  S(m),  (33) 

k= 0 
M-l 

J>(*  -  m)M*)  +  M*  -  mJftoW  =  S(m),  (34) 

k  —  0 
M-l 

Y,  9'(k  -  m)he(k)  +  g0(k  -  m)h0(k)  =  0.  (35) 

k=0 

Next  we  take  M-length  DFTs  of  (33)— (35)  to  obtain,  after 
changing  — u  to  u, 

l»e(«)|2  +  |g0(«)|2  =  1,  (36) 

IM«)|2+1M«0I2  =  1.  (37) 

Se(u)M«)  +  <7o(u)M«)  =  0.  (38) 

We  also  claim  that 


where  g.  h,  and  x  are  .V-periodic  sequences.  The  purpose 
of  this  subsection  is  to  exploit  the  formulas  (26)  and  (27) 
to  show  that  if  (4)  holds  and  if  h  is  given  by  (46)  below, 
then  (3),  (5),  and  (6)  also  hold.  The  derivation  is  adapted 
from  [2.  Section  3. A  ] 

Clearly,  G  and  H  can  be  identified  with  the  M  x  N 
matrices 


<j'(u)g0{u)  +  h'(u)'h0(u)  =  0.  (39) 

To  see  this,  consider  the  left-most  column  of  (3).  We  find 
that 

M-l 

Y  g(n  -  2m)g(—2m)  +  h(n  -  2m)/»(-2m)  =  6(n), 


'  <7(0) 

<7(1) 

g(N  -  1)  ' 

<7(  — 2) 

<7(-l) 

g{N  -  3) 

.  <7(2  -  A') 

<7(3  -W) 

<7(1)  . 

■  A(0) 

Ml) 

h(N  -  1)  ' 

M-2) 

M-D 

h(N  -3) 

h{2  -  N) 

h(  3-W) 

Ml)  . 

Regarding  (3)— (6)  as  matrix  equations  involving  (28)  and 
(29)  we  proceed  as  follows.  First,  consider  the  left-most 
columns  of  ( 4 )— (6) .  We  find  that 


m=:0 

(40) 

where  here  6  denotes  the  W-periodic  Kronecker  delta.  For 
n  =  2k  -t-  1  and  k  =  0, . . . ,  M  —  1,  this  becomes 

M-l 

y  <jo(k  -  m)ge(-m)  +  h0(k  -  m)h'(-m)  =  0. 

m=0 

Taking  M -length  DFTs  and  changing  —  u  to  u  yields  (39). 

We  now  claim  that  (36)-(39)  are  equivalent  to  It 

is  not  hard  to  see  that  (36)-(38)  are  respectively  equivalent 
to  (4)-(6).  However,  (39)  alone  does  not  imply  (3).  The 
reason  for  this,  as  will  become  clear  below,  is  that  we  must 
consider  both  the  even  and  the  odd  columns  of  (3). 


.v-i 


y  g(n  -  2m)g(n) 

=  Mm), 

(30) 

=  0 

.v-l 

y  M«  -  2m)h(n) 

n  =  0 

=  S(m), 

(31) 

Lemma  6:  Equations  (36)-(39)  together  imply  (3). 

Proof:  We  begin  with  the  following  observations 

First,  it  is  not  hard  to  show  that  (36)-(39)  imply 

IM«)I  =  l?o(ti)|  and  |M«i)|  =  jy«(u)|.  (41) 


Substituting  (41)  into  (36)  and  (37)  gives 


IM«)la  +  IM«)la  =  i. 

(42) 

IMM2  +  IMM2  =  l. 

(43) 

(Conversely,  it  is  also  true  that  (42),  (43),  (38),  and  (39) 
imply  (41),  and  hence  (36)  and  (37)  as  well.)  We  now  ob¬ 
serve  that  (42)  also  follows  directly  from  (40)  by  letting 
n  =  2k  and  taking  A/-length  DFTs.  To  obtain  (43)  di¬ 
rectly,  consider  the  second  column  of  (3)  (column  index  1). 
We  obtain 

V/-1 

r  g(n  -  2m)g(  1  —  2 m)  +  h(n  -  2m)h(  1  —  2m)  =  6(n—  1). 

m  =  0 

(44) 

Letting  n  =  2k  +  1  and  taking  A/-length  DFTs  yields  (43). 
We  ran  now  see  that  (36)-(39)  together  imply  (41).  which 
then  gives  us  (42)  and  (13)  Then  from  (39)  we  recover  (40) 
with  n  —  2k  +  1 ;  from  (42)  we  recover  (40)  with  n  =  2k, 
and  from  (43)  we  recover  (44)  with  n  =  2k  +  1.  Note  that 
(44)  with  n  =  2k  is  equivalent  to  (40)  with  n  =  2k  +  l. 
which  we  have  already  recovered.  □ 

Remark:  Similar  arguments  also  show  that  (3)  and  (6) 
alone  are  equivalent  (42),  (43),  (38).  and  (39).  This  was 
the  approach  used  in  (2].  With  either  approach  we  obtain 
(41).  on  which  we  focus  next. 

Equation  (41)  leads  us  to  impose  the  constraints 


Mu)  =  Mu)  and  Mu)  =  ~MM  (45) 
Lemma  7:  The  condition  (45)  is  equivalent  to 

h(n)  =  ( —  1  )n</(  1  —  n).  (46) 

Proof:  To  show  that  (45)  implies  (46),  substitute  (45) 
into  the  easily  verified  formula. 

Mu)  =  Mu)  +  /i0(u)e"-'#u,  (47) 

where  h  is  the  iV-length  DFT  of  h.  We  obtain 

Mu)  =  jo(u)  -  Mu)«~;^“, 

from  which  (46)  follows.  The  fact  that  (46)  implies  (45) 
follows  by  separately  considering  (46)  with  n  replaced  by 
2 n  +  1  and  by  2n  □ 

Theorem  8:  If  G  and  H  are  of  the  form  (26)  and  (27), 
respectively,  then  (30)  and  (46)  imply  that  (3)-(6)  hold. 

Proof:  Since  (46)  implies  (45),  we  clearly  have  (38) 
and  (39)  Next,  recall  that  (30)  is  equivalent  to  (36),  and  if 
(36)  holds,  then  (45)  implies  (37).  Finally,  since  ( 36 )-( 39) 
are  equivalent  to  (3)  (6).  the  theorem  follows.  □ 


Remark:  Solutions  of  (30)  are  easy  to  construct  because 
(30)  is  equivalent  to  (36).  Suppose  that  Mu)  's  arbitrarily 
given.  Without  loss  of  generality,  we  may  assume  that 
Mu)  has  been  scaled  so  that  |Mu)|2  <  1  for  all  u.  For  each 
u,  let  jo(u)  be  any  complex  number  such  that  |j0(u)|2  = 

1  —  |Mu)|2.  Then  (36)  holds  for  all  u.  If  we  let  h(n) 
be  given  by  (46),  then  G  and  H  will  satisfy  (3)-(6).  In 
the  next  section,  we  consider  the  more  delicate  task  of 
constructing  g(n)  so  that  in  addition  to  (30),  we  also  have 
^(n)  =  0  for  most  values  of  n. 

B.  Selection  of  Gq 

As  before  we  drop  the  subscript  0.  The  purpose  of  this 
subsection  is  to  find  a  solution  of  (30)  such  that  g{n)  —  0 
for  most  values  of  n. 

Let  A(m)  denote  the  sum  in  (30).  Clearly,  A  has  period 

A/.  Now  observe  that  as  m  takes  the  values  1, . \f / 2—  1 . 

the  number  M  —  m  takes  the  values  M  —  1 . V//2  +  1. 

Thus,  A(m)  =  0  for  m  =  1 . M/2  -  1  if  and  only  if 

A(m)  =  0  for  m  =  M/2  +  1  ,...,M  —  1.  Therefore,  (30) 
is  equivalent  to  saying  that  A(m)  =  0  for  m  =  0, . . . ,  M/2. 
We  now  impose  the  constraint  g(n)  =  0  for  p  <  n  <  N  —  1. 
where  p  <  M  and  p  is  even.  It  then  follows  that  (30)  can 
be  replaced  by 

p-i 

^MM<7(n  +  2m)  =  6(m),  m  =  0, . . . ,  M/2.  (48) 

n=0 

Now.  in  (48)  n  +  2m  <  (p  -  1)  +  2(M/2)  =  p  -  1  +  M  < 
<V  —  1,  since  we  assumed  p  <  M.  So,  the  sum  in  (48)  will 
automatically  be  zero  if  2m  >  p.  Thus,  we  may  rewrite 

(48)  as 

p- 1 

M n)<7(n  +  2m)  =  8(m),  m  =  0 . p/2-1.  (49) 

n=0 

Finally,  observe  that  n  +  2m  <  p  —  1  if  and  only  if  n  < 
p  -  2m  -  1 ,  and  since  m  <  p/2  -  1 ,  p  -  2m  -  1  >  1 .  Hence, 

(49)  becomes 

p-2m-l 

H  <?(")$(”  + 2m)  =  6(m),  m  =  0, . . . ,  p/2  -  1. 

n=0 

(50) 

Remark:  Obviously,  (50)  does  not  depend  on  N.  Now 
suppose  p  is  even  and  that  we  have  a  solution  to  (50) 
Suppose  also  that  N  and  J  are  such  that  2  <  p  <  N/2J . 
Clearly,  if  g,  is  given  by  (19),  then  (30)  will  hold  for  gr  If 
we  now  define  hj  by  (20),  Theorem  8  shows  that  (3)-(6) 
will  hold  for  Gj  and  Hj  if  they  are  defined  by  (14)  and 
(15),  respectively.  As  argued  in  Section  III- B,  we  obtain  a 
fast  wavelet  transform. 

We  now  observe  that  (50)  gives  us  only  p/2  nonlinear 
equations  in  p  unknowns.  Therefore,  we  must  impose  ad- 


ditional  consistent  constraints  on  $(0), . .  ,,g(p  -  1)  in  or¬ 
der  to  fix  a  solution.  The  linear  constraints  suggested  by 
Daubechies  [2j  are 

=°  for*  =  ° . P/ 2-1.  (51) 

'n=0  z  =  -l 

Theorem  9  ( Daubechies ):  The  constraints  (51)  are  con¬ 
sistent  with  (50).  In  fact,  we  can  require  that  g(n)  be  real 
valued. 

The  proof  of  this  result,  which  is  effectively  contained  in 
[2,  Section  4.B.],  is  sketched  in  the  Appendix. 

Remark  The  theorem  statement  only  guarantees  exis¬ 
tence  of  a  solution  to  (50)— (51 ).  In  fact,  the  proof  is  con¬ 
structive,  obtaining  solutions  by  factoring  certain  polyno¬ 
mials;  however,  this  is  not  the  only  way  to  obtain  a  solution 
in  a  given  instance.  One  can  attempt  a  direct  solution  by- 
hand  if  the  system  of  equations  is  small  (e  g  .  p  =  4)  or  by 
numerical  techniques  such  as  Newton's  method  otherwise. 
However,  see  also  the  remark  at  the  end  of  the  Appendix. 
Note  also  that  if  g(n)  is  real,  then  so  is  h{n).  Hence,  D 
can  also  be  viewed  as  a  mapping  from  IRV  to  IR,V 

C  Interpretation  of  (51) 

We  now  show  that  taking  k  —  0  in  (51)  suggests  that 
G  can  be  viewed  as  a  low-pass  filter  and  that  H  can  be 
viewed  as  a  high-pass  filter.  We  assume  that  (3)-(6)  hold 
with  G  and  H  given  by  ( 26)  ( 27 ) . 

Recalling  that  g(n)  =  0  for  n  =  p . N  -  1,  and  taking 

k  —  0  in  (51 )  yields 

.v- 1 

X>(n)(-U"  =  0,  (52) 

n  =0 

or.  equivalently, 

Af  —  1  Vf  -  1 

X  ?«(m)  =  X  $°(m) 

m  =  0  m  =  0 

In  other  words, 

iV-t 

X  9(n)  =  c.  (53) 

n  =0 

where  C  =  2  gf(m)  If  we  set  r(n)  =  1,  then  (53) 

implies  that 

v-i 

5^  g(ri  -  2m)x(n )  —  C  x(2m)  (54) 

n  =0 

In  other  words,  the  operator  G  simply  scales  and  subsam¬ 
ples  the  il.c.  signal  x(n)  Hence,  we  view  C  as  a  iow-pass 


filter.  Now  suppose  that  we  substitute  (46)  into  (52).  We 
find  that 

,v-i 

X  = 

n  =  0 

from  which  it  follows  that  if  x(n)  =  1,  then 
.v_i 

^  h(n  -  2m)x(n)  -  0. 

n=  0 

Since  H  blocks  the  d.c.  signal  .r(n),  we  view  H  as  a  high- 
pass  filter. 

We  now  show  that  the  constant  C  in  (54)  satisfies  C2  = 
2.  For  x(n)  =  1,  we  just  showed  that  Hx  =  0.  Thus,  (3) 
implies  j|G’Gx||2  =  ||x||2.  Next,  (4)  implies  ||G*G.r!j2  = 
||Gx||2.  Finally,  ||ij|2  =  S ,  and  by  (54),  ||Gx||2  =  C2.V/2. 

For  k  >  1,  (51)  simply  specifies  how  G  and  H  process 
certain  higher- variation  signals. 

Remark:  Daubechies'  original  reason  for  introduc¬ 

ing  (51)  was  related  to  regularizing  the  behavior  of 
the  continuous-time  analog  of  (18)  for  large  J  (2,  Sec¬ 
tion  3.B.], 

Appendix 

Proof  of  Theorem  9 

The  proof  is  self-contained  except  for  references  to  [2. 
Lemmas  4.2  and  4.4] .  Since  these  lemmas  are  statements 
about  certain  polynomials,  their  proofs  in  [2]  can  be  read 
independently  of  the  rest  of  [2], 

With  slight  abuse  of  notation,  let  G(;)  denote  the  poly¬ 
nomial  appearing  in  (51).  Then  (51)  is  equivalent  to  saying 
that  (1  +  ;)p/2  is  a  factor  of  G(z),  and  hence, 

G(z)  =  [i(l  +  -)]p/2F(x)  (56) 

for  some  polynomial  F  of  degree  p/2  —  1.  The  factor  of 
1/2  is  included  for  convenience  in  the  derivation  below 
Observe  that  the  coefficients  of  G  are  all  real  if  and  only 
if  the  coefficients  of  F  are  all  real.  We  now  require  that 
g(n)  be  real.  Thus,  it  suffices  to  prove  the  existence  of  a 
polynomial  F  of  degree  p/2  —  1  with  real  coefficients  such 
that  if  G(z)  is  defined  by  (56),  then  (50)  holds. 

To  begin,  first  recall  that  (50)  was  derived  from  (30), 
and  (30)  is  equivalent  to  (36).  If  we  let  g(u)  denote  the 
/V-length  DFT  of  y(n)  and  use  the  analog  of  (47),  we  find 
that  (36)  is  equivalent  to 

|<7(u)|2  +  |<7(u  +  A/2)|2  =  2.  (57) 

Now  observe  that  g(u)  —  G(e ;f)  when  £  =  —2 tu/,V  and  u 
is  an  integer.  If  (56)  holds,  then 

|G(e^)|2  =  [cos2(^/2)]p/2|F(ejf  )|2. 

a? 


Next,  if  F(z)  has  real  coefficients,  then  |F(eJ*)|2  also 
has  real  coefficients  and  is  real-valued.  Hence,  we  can 
write  |F(eJ<)|2  =  F(ej()F(e~i()  as  a  polynomial  of  de¬ 
gree  p/2  -  1  in  cos£,  or  equivalently,  in  sin2(£/2);  i.e..  we 
may  write  |F(eJ*)|2  =  P(sin2(£/2)),  where  P  is  a  polyno¬ 
mial  of  degree  p/2  -  1.  So,  if  (56)  holds  for  some  F  with 
real  coefficients,  and  if  (57)  holds  (which  is  what  we  are 
trying  to  prove),  then 

y"2P(l-y)  +  (l-y)"2P(y)  =  2.  (58) 

where  y  =  cos2(£/2),  £  =  -2iru/,V,  and  u  is  an  integer. 
Conversely,  in  [2,  p.  975,  Lemma  4.4]  it  is  proved  that  the 
polynomial 


Ply) 


P/ 2-1 

k=Q 


p/2-l  +  f\  fc 

k  'y 


satisfies  (58).  Furthermore,  P(y)  is  nonnegative  for  0  < 
y  <  1,  and  by  [2.  p.  972,  Lemma  4.2],  there  exists  a 
(nonunique)  polynomial  Q  of  degree  p/2  —  1  with  real  co¬ 
efficients  that  satisfies  |(?(ej{)|2  =  P(sin2(£/2)).  Thus,  if 
we  take  the  Q  of  Daubechies  and  obtain  the  g(n)  from 

G(z)  =  [*(1  +  .’))p/2QU), 

and  set  £  =  -2tru/,V  and  y  —  cos2(£/2),  then 

|y(u)|2  +  |y(u  +  .V/2)|2 

=  |G(e2<)|2  +  |C(e2«-'>)|2 

=  yW2p(l-y)  +  (l-y)p/2P(y) 

=  2. 

□ 

Remark:  In  general  there  are  p/2  choices  for  Q(z)  (and 
hence  (7(-))  The  coefficients  g(n)  (called  h(n)  in  [2])  given 
in  [2.  Table  1.  p.  980]  correspond  to  taking  <3(*)  to  have 
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Fig.  7.  Graph  of  x0(n)  defined  in  (22). 


Fig.  8.  Graph  of  x0(n)  defined  in  (23). 
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Fig.  9.  DWT  of  Fig.  7. 


Fig.  10.  DWT  of  Fig.  8. 
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Fig.  16.  Reconstruction  after  compression  of  signal  in  Fig.  14. 
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